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Abstract 


This  study  Addresses  the  maximization  of  the  circular  orbit  radius  of  an  earth-orbiting 
solar  sail  under  free  coning  motion.  The  objective  is  to  find  the  optimal  sail  settings 
producing  the  most  change  in  the  semi- major  axis  per  orbit.  Angular  orientations  with 
respect  to  sail  nutation,  precession,  mean  motion,  and  its  angular  momentum  control 
the  magnitude  of  the  solar  thrust  along  the  sailcrafl's  velocity  direction.  A  numerical 
search  scheme  uses  a  modified  Nevton-Rapson  iteration  method  to  identify  sets  of 
control  parameters  meeting  certain  optimality  conditions  that  produce  a  stationary 
value  in  a  selected  performance  index.  Such'scheme  displays  its  vulnerability  to  a  lack 
of  a  good  initial  guess.  Three  dimensional  perspectives  of  the  small  perturbation 
equation  describing  the  behavior  of  the  change  in  semi-major  axis  facilitates  the 
understanding  of  its  cyclic  nature  and  provides  an  excellent  tool  for  identifying  the 
various  locations  of  possible  maxima  (and  minima)  as  veil  as  slope-critical  area.  The 
perspectives  improves  the  initial  guess  for  implementing  maxima  search  schemes.  A 
test  case  demonstrates  the  location  of  particular  points  of  interest  vith  fev  search 
iterations.  /■  .■  -  '  0 
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u  Soi^f  Saiiifli 

The  Concept.  Solar  Sailing  is  a  method  of  propelling  an  object  through  space 
by  the  use  of  solar  radiation  pressure.  Issac  Newton  formulated  the  equations  that 
makes  this  method  realizable  and  most  feasible.  Pressure  acting  on  an  area  equates  to 
some  resultant  force  (force  -  pressure  z  area)  which.  if  not  acted  against,  would 
produce  some  type  of  motion  in  the  direction  of  that  force  (force  -  z  acceleration). 
One  can  think  of  solar  sailing  as  "sailing  with  the  sun".  Studies  accomplished  in  the 
past  have  shown  it  to  be  a  most  economical  method  of  interplanetary  travel.  The  key  to 
this  economy  is  in  the  nature  of  its  propulsion  system. 

The  Propulsion  System.  The  basic  propulsion  element  in  a  solar  sail  system  is 
a  highly  specularly  reflective  mirror-like  surface  that  creates  the  thrust  by  reflection 
of  sunlight.  The  physics  of  converting  photon  energy  to  spacecraft  motion  is  simply 
this:  the  solar  radiation  pressure  results  from  changes  in  the  momentum  of  incident 
photons  on  the  sails  reflective  surface.  The  higher  the  reflectivity  (i.e..  lower 
absorbtivity)  of  the  sail,  the  greater  would  be  this  momentum  transfer.  The  pressure 
on  the  ezposed  sail  surface  area  constitutes  the  "thrust"  on  the  vehicle  which  produces 
the  motion  (see  Fig.  1.1).  This  force  acts  to  increase  or  decrease  the  sail  total  energy. 
This  change  in  energy  allows  the  use  of  solar  sailing  as  a  means  of  space  travel. 
Controlling  the  angle  of  incidence  of  the  incoming  solar  raditaion  (photons)  amounts 
to  controlling  the  direction  of  this  resultant  force. 


The  Controls  System.  The  only  control  parameter  in  a  solar  sail  propulsion 
system  (with  constant  sail  surface  area)  is  the  "setting  angle"  of  the  sail  i.e.,  the  angle 
between  a  unit  vector  normal  to  the  shady  side  of  the  sail  and  a  unit  vector  in  the 
direction  of  the  sun  as  shown  in  Fig.  1.1.  The  effect  of  a  change  in  this  setting  angle  is 
a  change  in  the  direction  and  magnitude  of  the  resultant  solar  thrust  (force).  Since 
the  sail  changes  its  orientation  with  respect  to  an  inertial  frame  from  which  the 
orbital  parameters  are  referenced,  this  control  or  setting  angle  is  usually  described  by 
a  set  of  orientation  parameters. 

1.2  Maximising  the  Sami-Mninr  Axis 

Total  Energy  Approach.  In  order  to  increase  the  semi-major  axis  a,  the  total 
energy  (kinetic  plus  potential)  of  the  sailcrafl  must  be  increased  from  the  energy  of 
the  original  orbit  to  some  higher  energy  level.  The  rate  of  energy  increase  is  equal  to 
the  rate  that  work  is  done  on  the  sailcrafl  by  the  thrust  produced  by  its  propulsion 
system  (solar  sails  themselves).  It  may  seem  advantageous  and  optimal  to  direct  the 
thrust  such  that  the  rate  of  work  done  on  the  sailcrafl  is  always  at  or  near  a  maximum. 


This  is  not  necessarily  true.  For  circular  orbits,  maximizing  the  semi-major  axis  is 
essentially  maximizing  the  orbital  radius.  For  small  perturbations,  this  assumption  is 
valid.  Maximizing  the  rate  of  energy  increase  is  realized  for  this  sailcrafl  Then  the 
thrust  component  in  the  velocity  vector  is  maximized.  Analyses  have  shown  that  the 
orientation  of  the  sailcrafl' s  orbit  with  respect  to  the  earth-sun  line  determines  the 
maximum  posible  rate  of  this  enrgy  increase.  [  Ref:  4]  All  strategies  to  increase  the 
semi-major  axis  are  indeed  energy  increase  strategies. 

1.3  Problem  Statement 

A  Solar  sail  being  in  a  photon  rich  environment,  inevitably  experience  a  perturbational 
force  which  affect  the  orbital  characteristics.  These  effectss  are  indeed  dependent 
upon  the  sails  orientation  with  respect  to  the  solar  radiation  source.  To  effectively  use 
the  resulting  thrust  created  by  this  radiation  pressure  to  provide  changes  in  the  orbital 
parameters,  the  orientation  must  be  known.  Random  orientation  will  result  in  random 
orbital  parameter  changes.  For  space  travel,  interest  lies  in  the  maximum  changes  in 
orbital  parameter  the  semi-major  axis.  Certain  combination  of  solar  sail  setting  angles 
can  provide  the  best  change  for  a  given  cost  parameter  ....  be  it  travel  time  or  the 
number  of  revolutions  to  acquire  a  desired  change. 

14  Objective  and  Scone. 

The  objective  of  this  study  is  to  determine  what  set  of  control  angles  would  provide  the 
optimal  change  in  the  semi-major  axis  of  an  earth-orbiting  solar  sailcralt.  The  scope 
of  this  effort  is  to  investigate  the  nature  the  perturbation  of  the  semi- major  axis  and 
determine,  if  any.  the  set  of  control  parameters  that  produce  the  maximum  change  in 
this  orbital  parameter  in  one  revolution. 


Technique  This  effort  incorporates  the  deserption  of  semi-major  axis 
perturbation  function  and  the  determination  of  conditions  that  it.  The 

technique  of  dynamical  systems  optimization  is  employed  as  presented  by  Bryson  and 
Hb  [Ref:  2d) 

Tools.  This  approach  requires  the  identification  of  the  perturbation  function 
and  its  corresponding  surface  followed  by  the  construction  of  a  computer  program  to 
evaluate  the  sensitivity  of  the  orbital  parameters  to  small  perturbational  forces  and  to 
search  for  a  set  of  control  angles  satisfying  a  given  performance  index.  This  index  is 
established  as  the  "maximum  change  in  semi-major  axis  per  orbit*. 

1.6  Summary  of  Report. 

Section  II  provides  some  background  information  into  the  previous  studies  on 
solar  sailing  and  related  topics.  Their  results  are  summarized.  A  description  of  the  solar 
sail  model  used  in  the  study  along  with  the  basic  and  simplifying  assumptions  that  led 
to  development  of  the  solar  sail  coning  phenomenon  follows  A  short  rationale  for  the 
neglect  of  any  shadowing  effects  is  presented.  The  general  perturbation  equation 
derived  by  an  earlier  researcher  is  examined  for  a  special  resonance  case  and  is 
evaluated  for  three  distinct  sailcrafl  motions,  spinning,  coning,  and  tumbling. 

Section  III  discusses  the  use  of  surface  representations  of  the  perturbation 
function  Aa  as  a  method  to  identify  the  regions  of  wwiM  and  ■»««««»*  The  dynamic 
system  optimization  that  follows  it  bases  its  region-of-search  on  these  surface 
representations.  The  approach  to  the  solution  via  the  gradient  method  is  presented. 

Section  IV  provides  an  insight  into  the  nature  of  the  perturbation  equation 
by  graphically  portraying  the  surface  of  the  perturbation  function.  These  surfaces 
represent  the  fractions  behavior  from  a  three-dimensional  perspective  at  critical 
setting  angles  and  provides  the  best  initial  guess  for  the  search  scheme 


Section V  presents  the  conclusions  drawn  from  the  entire  effort.  It  includes 
the  lessons  learned  and  suggestions  for  further  study. 

The  Appendices  provide  a  detailed  look  at  the  eclipsing  (shadowing) 
phenomenon(Appendix  A)  along  with  a  rationale  for  employing  the  numerical 
technique  in  lieu  of  the  analytical  (Appendix  B).  The  Fortran  Programs  used  to 
generate  the  numerous  surfaces  and  to  determine  the  optimum  control  settings  are 
included  in  Appendix  C  and  D,  respectively.  Appendix  E  provides  a  sample  output  of  a 
test  case  to  show  the  search  mechanism  using  the  Stirling  approximation  technique. 


II.  Background 


2.1  Previous  Efforts 

There  hove  been  severe!  studies  in  the  pest  that  have  placed  interest  in  the  solar 
sailing  concepts  and  in  the  solar  radiation  induced  orbital  perturbation  of  space 
structures.  Studies  in  the  solar  sailing  concept  tend  to  deal  vith  space  travel  and  the 
determination  of  control  lavs  to  provide  optimum  changes  in  orbital  parameters  such 
as  semi-major  axis  and  inclination  to  attain  escape  from  a  planet's  gravitational  field. 
Studies  in  the  orbital  perturbation  due  to  solar  radiation  have  generally  employed 
numerical  and  expansions  techniques  to  determine  orbital  motion  due  to  the 
"formidability"  of  multi-first  order,  nonlinear  coupled  differential  equations. 

Tsander.  [Ref:  18]  .  is  considered  the  "father"  of  solar  sailing.  His  serious 
investigations  of  the  solar  sail  problem  for  spaceflight  demonstrated,  in  principle,  the 
feasibility  of  making  interplanetary  flights  with  the  aid  of  solar  pressure. 


solar  sail  (furled) 

i 


solar  sail  (uafarled) 


[Ref:  31 


Fig  2 Jk:  Furling/UnAirling  Method 


(direction  of  motion) 


(vith  sail  set  for  approaching  sun) 


Fig.  2. IB:  Optimum  Angle  Orientation 


(Ref:  5) 


Garvin.  (Ref  3].  vhose  vork  in  establishing  the  concept  of  solar  sailing  as  a 
practical  means  for  space  vehicle  propulsion,  considers  the  "furling  and  unfurling 
technique  to  increase  the  altitude  of  a  sailcrafl  and  escape  the  earth's  gravitational 
field.  Figure  2.1A  sheers  this  method.  The  objective  is  to  a  as  i  mire  change  in 
semi-major  axis  as  the  sailcrafl  travels  avay  from  the  sun  and  to  minimize  the  change 
as  it  approaches  the  sun  is  the  apparent  rationale  for  this  approach.  This  study  states 

that  an  optimum  tilt  angle  is  8  »33*  for  maximizing  the  semi-major  axis.  From  figure 
2. IB  this  angle  vill  result  in  the  largest  thrust  along  the  velocity  vector  (i.e..  along  the 
direction  of  motion ). 

Fimnle  (Ref:  41 .  analyzes  the  use  of  solar  sails  for  an  earth  escape  trajectory 
from  a  circular  orbit  using  the  maximum  time  rating  energy  increase  approach. 
Trajectory  analyses  for  two  different  orbit  geometries  (see  Figure  22)  reveals  that 
orbits  parallel  to  the  solar  radiation  (Orbit  A)  results  in  a  radical  orbit  eccentricity 
change.  Orbits  normal  to  the  solar  radiation  (Orbit  B)  result  in  a  slovly  changing  orbit 
eccentricity  or  quasi-circular  trajectories.  Hence.  Orbit  B  enhances  the  change  in 
orbit  semi-major  axis:  launching  into  this  orbit  is .  hovever.  a  difficulty 


Fig  2  2:  Orbit  A  Geometry 


fltef:  4) 


(Ref:  4) 


Fig  22B:  Orbit  B  Geometry 
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Sands.  [Ref  .  151  .  studied  the  similar  use  of  solar  sails  for  escaping  the 
gravitational  influence  of  a  given  planet  employing  the  total  energy  change 
approach.  Results  shov  that  placing  the  sailcraft  in  an  initially  circular  orbit 

spinning  at  a  rate  (J  -  1/2  period  constitutes  an  inefficent  escape  maneuver.  An 
alternative  suggested  was  an  eliptical  orbit  vith  energy  level  as  the  original  starting 
orbit. 


Cotter.  [Ref:  27] ,  introduces  a  very  important  number  called  the  "lightness"  of 
the  sail.  This  lightness  is  the  ratio  of  the  maximum  solar  radiation  pressure  force  to  the 
solar  gravity  force  on  a  given  sailcraft.  Since  these  tvo  forces  are  functions  of  the 
inverse  square  lav  in  vhich  its  value  varies  inversely  as  the  square  of  the  distance 
from  the  sun.  this  ratio  independent  of  this  distance  from  the  sun.  This  maJtes  this 
ratio  uniquely  a  directly  proportional  measure  of  the  inertia  of  the  sailcraft  (i.e..  the 
larger  the  ratio,  the  greater  the  acceleration  of  the  sailcraft). 

This  study  discusses  the  feasibility  of  solar  sailing  and  demonstrates  it  application  for 
interplanetary  travel.  Travel  limitations  are  in  time  and  temperature.  The 
temperature  limitation  (operating  temperature  envelope  of  the  sailcraft)  dictates  a 
prudent  use  of  near-sun  trajectories;  the  time  limitations  restrict  planetary  travel  to 
the  nearer  portions  of  the  solar  system.  The  advantage  of  having  unlimited  "fuel" 
supply  makes  this  solar  sailing  concept  most  promising 
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London.  [Ref:  12]  .  provides  some  insight  into  the  motion  of  a  solar  sailcraft 
with  constant  sail  settings.  The  study  employs  the  logarithmic  spiral  trajectories  as  a 
solution  and  concludes  that  it  (solution )  is  not  optimum.  It  encourages  the  use  of  other 
trajectories  for  more  efficient  sail  utilization. 

Tsu.  [Ref:  19] .  evaluates  the  neccessity  and  importance  of  low  mass  to  area 
ratios  for  solar  sail  design  along  with  the  travel  time  and  optimization  of  the  sail  tilt 

angle  0  for  interplanetary  trips  with  minimum  time  as  a  constraint.  The  optimum  sail 

setting  for  such  application  is  found  to  be  dependent  upon  the  acceleration  <X  of  the 
sailcraft  (due  to  radiation  pressure  acting  on  the  sail  area)  and  the  sun's  gavitational 
acceleration  at  earth's  orbit. 

The  time  of  travel  between  radii  r0  and  rf  is 

t  =  1/3  3/2  ~  rf 3/2 )  e  (aft/«  -  cos38)  1/2  (2.1.a) 

r0  <x,/2  si&0  cos29 

where  r0  -  initial  orbit  radius, 
rf  -  final  orbit  radius. 

Sq  -  sun's  gravitational  acceleration, 

a  -  sailcraft  acceleration  due  to  solar  pressure. 

0  -  angle  of  incidence. 

For  shortest  time. 

(2.1. b) 


from  which  the  optimum  tilt  angle  0  0pt  can  be  found.  Tsu  plots  this  out  and  finds  that 

for  Eq*  0.  the  optimum  tilt  angle  0Qpt  is  approximately  33°  This  is  the  same  result  found 
by  Garwin  (Ref:  31. 


d  (a^/cc  -  cos  3  8) 1/2 
dO  sin9  cos20 


Van  der  Ha.  and  Modi.  (Ref:  21A1 ,  have  done  the  most  extensive  analytical  studies 
of  orbital  perturbation  caused  by  solar  radiation.  This  study  considers  a  more  realistic 
solar  radiation  model  (than  the  model  presented  in  section  2.2)  and  analyze  three  special 
orientation  cases  vith  a)  fixed  angle  with  respect  (vrt)  to  the  solar  radiation  source,  b) 
fixed  angle  vrt  an  inertial  freme  and  c)  a  general  fixed  orientation  vrt  the  earth. 


Fig.  2.3:  Switching  Points  Configuration 

[Ref:  21A] 

In  another  study  (Ref:  21C1.  these  same  authors  employed  on/off  switching  strategies  to 
increase  the  semi-major  axis.  Such  strategies  involve  instantaneous  switching  controls 
which  in  effect  turns  the  solar  pressure  force  off  during  certain  parts  of  the  orbit.  The 
strategy  calls  for  "switch-on"  when  the  saiicraft  is  (roughly)  near  the  earth-sun  line 
(point  1)  and  "switch-off  at  (point  2)  as  shown  in  Figure  2.3  (circular  orbit  shown). 

Although  instantaneous  on/off  switching  may  be  technically  impractical,  it  is 
theoretically  effective  in  the  sense  that  the  rate  of  change  of  the  total  energy  is  always 
positive  during  the  on-phase  since  the  component  of  the  perturbing  force  along  the 
instantaneous  velocity  vector  is  positive. 


tonkins  [Ref:  101 .  is  the  only  individual  to  date  who  has  explored  the  behavior  of 
the  coning  solar  sail.  The  coning  motion  of  a  given  sailcraft  greatly  affects  the 
magnitude  and  direction  of  solar  radiation  pressure  force.  His  study  generated  the 

perturbation  algorithms  that  determines  the  changes  in  orbital  parameters  (a,  e,  i,  Q)  as 
a  function  of  specified  sail  setting  angles  which  relate  the  sailcraft  orientaition  with 
respect  to  an  inertial  frame. 


2.2  The  Coning  Solar  Sail. 


Solar  Sail  Model.  In  general,  a  solar  sailcrafl  must  meet  one  critical 
criterion  before  it  can  literally  sail"  anywhere.  The  mass-to-area  ratio  is  this  critical 
parameter  that  must  be  established.  It  is  the  ratio  of  the  total  solar  sail  area  to  the  total 
mass  of  the  sailcrafl  (sail  ♦  structure  ♦  electronics  *  payload).  Studies  have  shown  that 
the  mass/area  must  be  very  low.  This  is  quite  apparent  from  the  equation  for 
acceleration  due  to  solar  measure.  The  solar  acceleration  due  to  solar  pressure  is 
expressed  as 

a  solar  pressure  “  Po  A/M  ^.2) 

where  p0  -  solar  radiation  pressure  at  earth  s  orbit  radius. 

A  -  Area  of  Solar  sail  and  M  -  total  mass  of  solar  sailcrafl. 


m  m  '?>'  ,  wm 

. 

- .... - ///y/z/wga 


Fig  2.4:  Solar  Sail  Model 


(Ref:  10} 
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Typical  values  for  p0  is  about  9  x  10"*  dyne/cm2.  For  any  appreciable  amount  of 

acceleration  for  space  travel.  A/M  must  be  high  (or  its  reciprocal  M/A  be  very  low). 
This  essentially  requires  that  sail  area  A  be  much  larger  vith  respect  to  the  sailcraft 
mass.  [Ref:  19] 

The  solar  sail  model  used  in  this  effort  is  modeled  as  a  perfectly,  specularly  reflective 
(both  sides)  thin,  flat  rigid  plate  as  shown  above  in  Figure  2.4.  Recall,  a  highly 
specular  surface  means  that  all  incident  photons  are  reflected.  Hence,  the  solar  sail  is 
assumed  to  have  a  reflectivity  of  1.0  which  implies  that  all  incident  radiation  is 
converted  to  thrust.  Having  both  sides  with  the  same  sail  characteristics  means  that  all 
surfaces  exposed  to  direct  solar  radiation  will  be  contributing  to  the  overall  thrust  of 
the  vehicle.  The  A/M  ratio  for  this  model  is  assumed  large  enough  to  make  its  design 
feasible  for  solar  sailing. 

Basic  Assumptions.  The  basic  assumptions  made  in  the  course  of  this  effort  ana 
in  the  development  of  the  perturbation  solution  which  follows  are  enumerated  below 
categorically: 

A.  For  satellite  motion,  it  is  assumed  that . 

la.  the  solar  sailcraft  is  initially  injected  into  a  high  altitude  circular 
orbit; 

2a.  the  gravitational  field  is  central; 

3a.  perturbations  to  the  central  force  gravity  field  due  to  oblateness  and 
the  gravitational  effect  of  other  celestial  bodies  (e  g.,  moon,  sun) 
are  neglected;  also  neglected  are  any  magnetic  distorbances 
4a.  the  solar  sail  motion  is  torque-free;  therefore,  the  angular 

momentum  H  is  conserved. 

5a.  the  angular  velocity  vector  U  is  not  necessarily  colinear  with  the 


angular  momentum  vector  H; 

6a  for  eny  orientation,  dfl/dt  -  0  -»  nutation  angle  8  -  constant: 

7a.  for  any  orientation.  dV /dt  -  0  -»  precession  rate  V  -  constant: 

8a.  orbit  plane  is  considered  inertial: 

9a.  perturbations!  changes  are  small;  and 

10a.  the  body  spin  rate  0(  is  large  enough  for  spin  stabilization. 

B.  For  the  solar  sail  it  is  assumed  that.  .. 

lb.  the  spinning  solar  sail  is  an  axi-symmetric  rigid  body  vhose 
sail  surface  does  not  deform  under  loads: 

2b.  the  solar  sail  surface  characteristic  are  essentially  homogenous 
and  time  invariant; 

3b  the  solar  sail  is  perfectly  specular  reflective  on  both  sides;  and 

4b.  the  mass/area  ratio  «  1 .0. 

C.  For  solar  radiation,  it  is  assumed  that.  .. 

lc.  reflection  of  earth  surface  is  neglected;  the  sun  is  the  only  source 
of  radiation;  and 

2c.  no  shadovina  or  eclipsing  of  sailcraft  occurs.  This  is  due  to  the 
high  altitude  and  the  inclination  of  the  orbit  in  case. 


Torgue-Free  Motion  and  the  Coning  Phenomenon.  This  section  considers  the 
torque  motion  of  the  solar  model  vhich  is  axis-symmetric  and  has  a  principal  Moment 

of  Inertia  [A*  A'  C]  vith  A'  along  the  symmetric  body  axes  and  C  along  the  bj  axis 
perpendicular  to  the  plane  of  the  sail.  The  general  equation  of  torque-free  motion  is 

£  Mjjj-  *  dH/dt  *  0  (2.3) 


since  (by  definition)  Torque-Free  ^  moment  M  equals  zero.  R  must  be  constant  in  both 
direction  and  magnitude;  i.e.. 


H  - 


(2.4) 


▼here  Hq  is  the  initial  angular  momentum  of  the  sailcraft  about  its  center  of  mass 

(vhich  is  assumed  to  coincide  ▼ith  its  geometric  center).  The  direction  of  H  is  fixed  in 
inertial  space  and  can.  therefore,  can  be  arbitrarily  referenced  to  an  inertial  frame. 
Figures  2.5A  &  2.5B  depict  this  solar  sailcraft  ▼ith  respect  to  three  different  reference 
frames:  (Ref:  2) 


a)  Body  Fixed  Frame:  (bj  b2  b^) 

b)  Orbital  (inertial)  Frame:  (i  j  k) 

c)  Body-Centric  (inertial)  Frame:  (T  J  R) 


Using  the  same  angle  notations  as  found  in  [Ref:  10) .  Hfl  can  be  expressed  in  terms  of 
Euler  angles  as  folio  vs: 

Rq  =  H0sin8sinoct  bj  +  HQsin0  cosat  62  ♦  HgCOseb^  (23a) 


However .  for  principal  axis,  this  fiQ  equates  (by  definition)  to  the  following  expression: 
Rq  •  A'  1*1 1  bj  ♦  A'  UI2  t>2  *  C'  Uj  [Ref:  9)  (2.3b) 

where  { A'  A’  Cl  are  the  Principal  Moments  of  Inertia.  In  matrix  notation,  this  can  be 
expressed  readily  as  follows: 


» 

H0  sine  sinat 

r*li 

> 

£ 

9C 

O 

M 

H0  sin8  cosat 

i 

*2 

* 

A  “52 

H0cose 

. 

lb3. 

c  ^3 

(26a) 

(2.6b) 

(2.6c) 


Equating  term  per  term  and  solving  for  the  angular  velocity  about  each  body  axis . 


tJjj!  _  Hn  sine  sincct  _  d 

A  dt 


H0  sinBjinat 


(2.7a) 


uK,  _  Ha  sine  sinat  ufc,  _  d 

Hn  sine  cosat 

A' 

dt  l 

A 

l 

“b3  =  5oi225 

“»  “b3  =  i 

Ha  cos  8 

>0 

C’ 

dt 

C’ 

(2.7b) 


(2.7c) 


The  Euler  equations  for  dynamical  motion  (as  applied  to  this  sailcraft)  are  given  as 


follows:  [Ref:  11] . 

ZMbi  -  A  (Jjjj  -  (A  -  C )  (*>t>2ub3  •  0 

(2.8a) 

^*^b2  “  A  wb2  "  "  A  )  wb3wbl  “  ® 

(2.8b) 

ZMb3  -  C'tdb3  -  0 

(2.8c) 

Applying  the  respective  values  for  ths  tiaae  derivatives  of  the  angular  velocities,  the 
results  are  as  follow 

A  d  Hojin0_sinat  +  (C  -  A  )  H02  sin0  cos0  sinat  ,  0  (2.9a) 

dt  A  AC' 

A'  d  H0  sm|  cos  at (A~  -  C)  H02  sin0  cos0  sinat  „  q  (2  9b) 

dt  A  AC 

A'  d  Ho_COS0  m  q  (2.9c) 

dt  q- 

From  Equation  (2.9c).  the  derived  conclusion  is  that  Ho(cos0)/C’-  Constant;  hence. 

0  *  0O  *  constant  (2.10) 

since  H0  and  C'are  established  as  constants.  Hence  the  nutation  angle  or  the  coning 
angle  remains  fixed. 

Using  this  fact  in  Equation  (2.9b)  and  carrying  out  the  differentiation,  the  result  is  as 
follow:  The  magnitude  of  the  angular  velocity  is . 

a  -  (A’  -  C)  HOCOS0O  (2.11) 

AC 

from  which  is  known  that  (recall,  the  Moment  of  Inertia  terms  and  H0  are  constant).. 

a  s  the  body  angular  velocity  -  constant.  (2.12) 

The  precession  is  found  from  the  components  of  the  angular  velocity  about  body  axis 
b}.  this  happens  to  be  (as  reference  from  Figure  2.5A) 


Uhi 


a  ♦  v  cos0  =  cos0/C 


(213) 


Solving  for  V  (with  0(  as  expressed  above),  the  result  is 


V  -  EL  /  A'  -  constant 


Fig.  2.6:  Velocity  Vectors 

Solving  for  the  total  angular  velocity  U  is  the  vector  sum  of  a  and  V  (see  Figure  2.6). 

u  -  a  ♦  v  (2.15) 

and  would  be  at  a  set  angle  6  from  the  R  direction.  Vith  V  (precession  rate)  fixed 
along  R  direction  and  a  (body  spin)  fixed  along  body  b^  direction,  this  phenomenon 
displays  the  vector  U  as  "sweeping  out”  a  cone.  The  resultant  motion  of  this 

axis-symmetric  rigid  body  in  which  the  angular  velocity  vector  (D  is  not  colinear  with 
the  angular  momentum  vector  H  is  known  as  "coning".  This  behavior  is  shown  in 
Figure  2.7. 


[Ref:  2.10] 


Fig.  2.7:  Co&ing  Motion 


The  coning  angle  (as  referred  to  in  this  ten)  is  8.  It,  too.  displays  the  "sweeping"  cone 
phenomenon.  This  is  easily  seen  with  a  projection  of  the  body  b  j  axis  on  the  IJ  plane. 
This  is  accomplished  in  Figure  2.8  .  Line  OX  is  this  projection  and  it  precesses  about  the 
K  axis  at  a  rate  V.  the  precession  rate. 


The  Orientation  Angles  (8.  Tl.  c.  6-dr).  The  coning  angle  8  (as  discussed 


previously)  is  the  angle  between  the  angular  momentem  vector  H  and  the  Solar  Sail 
body  axis  (Ref:  Fig.  2.3B].  The  Angular  momentum  vector  H  is  referenced  to  the 
orbital  plane  via  angles  11  and  C  as  shown  in  Figure  2.9A  below: 


Fig.2.9A:  T|  and  £  Orientation 


[Ref:  10] 

The  IJK  Frame  is  defined  to  be  sail  body-centered  with  the  I  axis  in  a  plane 
perpendicular  to  the  orbital  plane.  This  study  refers  heavily  on  the  difference 

between  two  phase  angles,  $  and  This  difference,  suitably  called  PA.  equals  6- dr 


\|r  Phase  Angle 
_  -  B  Coning  Angle 

^  V  Precession  Rate 


Fig.  2.9C:  Angle  \|r  Orientation 


With,  8  -  T1  -  C  ■  0*.  the  angle  can  be  best  seen  for  understanding  as  is 

done  in  Figure  2.9D  below: 


Fig.  2.9D:  Angle  ($-^)  Orientation 
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2.3  Shadowing  Of  Eclipsing  Effect. 


Eclipsing  or  shadowing  of  the  solar  sail  was  not  considered  in  the  formulation  of  the 
perturbation  solutions.  Investigation  of  this  shadow  phenomenon  lead  to  the  conclusion 
that  for  a  given  orbit  inclination  and  altitude,  shadow  effects  can  be  ignored.  From  a 
general  geometrical  perspective,  one  can  safely  conclude  that  shadowing  will  occur  if 
the  orbit  lies  in  the  ecliptic  (earth-sun)  plane.  This  is  theortically  true.  The  interesting 
question  is  this:  At  what  orbit  inclination  and  altitude  is  shadowing  a  problem?  Stoddard 
(Ref:  25)  evaiutes  this  phenomenon  and  arrives  at  a  set  of  equations  that  determine  the... 

a)  minimum  inclination  before  shadowing  effects  can  be  considered  dominant 
and  b)  the  altitudes  a  satellite  must  be  at  to  avoid  shadowing. 

Filler  (Ref:  6)  circumvents  this  entire  issue  by  studying  solar  sailing  at  orbits 
perpendicular  to  the  ecliptic  plane.  This  approach  surely  removes  any  uncertainties  of 
shadow  effects  completely.  Solar  sails,  as  dependent  as  they  are  to  the  amount  of 
incident  radiation  for  propulsion,  can  be  so  placed  in  an  orientation  that  minimizes  this 
eclipsing  phenomenon. 

Figure  2.10  below  shows  the  shadow  and  no-shadow  situations.  From  this,  one  can  state 
that  there  does  exist  some  physical  limit  to  the  shadowing.  Stoddard  explains  and 
develops  a  method  to  analyze  this  situation.  He  presents  his  arguments  using  the 
cylindrical  shadow  theory  and  concludes  that  the  critical  inclination  for  a  given  orbit 
radius  is... 

is  cos_1lRe/Rsl  (2.17) 

Stoddard  claims  that  for 

(2.18a) 


Rs  <  R0  esc  i  =»  an  eclipse  will  take  place. 

R,  >  Re  esc  i  =»  an  eclipse  will  ggt  take  place. 


(2.18b) 


His  arguments  on  this  issue,  along  with  a  determination  for  the  duration  of  shadowing  (if 
any),  are  provided  in  Appendix  A.  The  results  are  interesting.  They  show  the 
comparison  of  a  system  with  and  without  shadow  effects  accounted  for.  The  hoti/im  tine  is 
this:  For  low-earth  orbit  solar  sails  within  the  shadow  limits,  shadowing  effects  are 
dominant  and  mum  be  considered.  By  increasing  the  orbit  inclination  and  the  attitude,  a 

satellite  can  circumvent  this  shadowing.  This  shadowing  phenomenon  and  its  effects 
are  circumvented  in  this  study  by  employing  non-ecliptic  planes  and  orbit  altitudes 
beyond  the  shadow  regime.  Vith  this  understood,  the  tasks  of  analyzing  solar  sail 
behavior  in  a  photon  rich  environment  is  pursued  here. 


2  4  General  Perturbation  Equation. 


Continuing  Work  This  study  continues  the  work  initiated  by  Jenkins  (Ref:  10) 
on  the  study  of  the  "Orbital  Motion  of  Coning  Solar  Sails"  in  which  he  succeeded  in 
developing  a  set  of  algorithms  for  determining  the  perturbations  in  certain  orbital 
parameters  (a.  e  1.  Q )  as  a  function  of  the  solar  sail  settings.  These  solar  settings  are 

lust  the  ortr  uoa  angles  (8,  7),  (,  $-40  of  the  solar  sail  with  respect  to  the  solar 

vector  5  Since  these  setting  angles  control  the  behavior  of  the  solar  sailcraft.  it  was  of 
particular  interest  to  determine  what  angle  (s)  is/are  dominant  in  producing  the 
greatest  changes  in  the  orbital  parameters)  selected.  This  set  of  algorithms  was 
derived  from  the  Lagrangian  Planetary  Equations  and  are  presented  below  in  a  more 
modular  fashion  than  found  in  [Ref:  10]. 

Perturbation  Equations. 


Aa  =  f  (da/dt)  dt 

Jo 

(2.19  a) 

Ae  =  {  (de/di)  dt 

Jo 

(2.19b) 

Ai  =  }T  (di/dt)  dt 

(2.19c) 

AQ  =  T  (dQ/dt)  dt 

(2. 19d) 

For  circular  orbits  (e  ■  0).  the  Langrangian  Planetary  Equations  reduce  to  the 
following: 


(2.20a) 


da  /dt  .  2P(b,.S)2(b,.9) 
n 

de/dt  =  D  (b,»5)2(b,aQ)  sin  f  +  2  D  [(E,e5)2(b,a9)  cos  fl  (220b) 

na 

di/dt  =  D  (b,*5)2  (b,*W)  cos  f  (2.20c) 

na 

dQ/dt  =  D  (b,*S)2  (b^W)  sin  f  (2.20d) 

na  sin  i 

vith  D  =  Sail  constant  -  3k  kg/m 
f  -  true  anomaly 

bg  -  body  vector  b-j  expressed  in  orbital  reference  frame 
•  s  dot  product 
§  -  solar  vector  (direction  of  sun) 

(GVW)  -  (Radial.  Tangential,  Normal)  components  of  thrust  per  unit  mass. 

One-To-One  Resonance  Perturbation.  In  the  context  of  this  study,  the 
one-to-one  resonance  occurs  vhen  the  mean  motion  of  the  satellite  about  its  planet 


equals  its  precession  rate  about  it  angular  momentum  vector  H.  If  the  integration  in 
Equation  2.19  is  taken  over  one  orbital  period  with  the  orbital  motion  (n)  equal  to  its 

precession  rate  (V),  the  following  expressions  result: 

1)  Aa  =  a  L5DTP/u  [da,  +  da2+  da3  +  da4  +  da5  +  da6]  (2.21) 


where 

da,  »  !4  D,2  I(p2+3p4)  cos($-y)  -  (3pl+p5)  sin($-^)  1  (2.22a) 

da2  =  '/iD,  D2  I(p2-p4)sin($-tif)  +  (p5-pi  )cos($-4r)l  (2.22b) 

da3  -  2  D,  D3  [  p6  cos($->|0  ♦  P3  sin($->|0  ]  (222c) 

da4  *  2  D2  D3  l  ^3  cos($-^)  +  p6  sin($-^)  1  (2.22d) 

da5  *  !d  D2  D2  [  (3P2+P4)  cos($-*)  -  (£1  +3p5)sin(<|>-*)  ]  (2.22e) 

da6  =  D3  2  [  (p2  +P4)  cos($-*)  -  (pi  ♦  p5)  sin($-*)  ]  (2.22f) 

2)  Ae  =  a^DTP/8ji  [de,  +  de2  +  de3  +  de4  +  de3]  (2.23) 


where 

de,  =  D,  2  [  -p3  sin(2($-*))  ♦  p6  (6  ♦  cos(2($-*)))  ]  (2.24a) 

de2  =  2  D,  D2  [  p3  cos(2($-4r))  ♦  P6  sin(2(<t>-^))  1 


(2.24b) 


de3  - 2 D,  D3  1(02+04) cos(2tt-*)M01+05) sin(2(^-^))  +  6^)1  (224c) 
de<  -  2  D2  D3  1(02+04)  sin(2(^-*))+(01+05)  cos(2(^))-605  I  <224d) 
de5»12D32p6  (224e) 


3)  A i  -  a4DTP/2ji  [  di,  ♦  di2  ] 


▼here 


(223) 


(226b) 


(2.27) 


di,  *  208  D3  ID,  oos(^-^)  ♦  D2  sin($-*)  ]  (226a) 

dij  «  07  l  ( D,  2+.75  (D2)2+  D3  2)  sin(^-x|r)  +.5  D,  Dacos(^Mr)  1  (226b) 

4)  AQ  -  a*  DTP  /(2ji  sin  i)  [dQ,  ♦  dQ2]  (2.27) 

▼here 

dQ,  =  0,1(  25  D,2  +  .75  D22  +  D32]  cos  (*  -  *)  -.5  D,  D2  sin  (♦  -  ♦)]  (228a) 

dQ2  -  20,  D3  (D,  sin  (4-^)  -  D2  cos($  -  4#)];  (228b) 


The  other  parameters  are  as  follow: 

D,  -  sin  0  cos  (  cos  i 

D2  »  sin  0  sin  C  cos  T)  cos  i  ♦  sin  0  sin  T)  sin  i 


(229a) 


(229b) 


’*  iA  .  ■  '  ■  . 


jmr.o.vo.omtr. 


D5  -  cos  6  sin  C  sin  T)  cos  i  ♦  cos  0  cos  T)  sin  i 


(229c) 


tad 


pi  -  sin0cos< 

(2.30a) 

02  -  sin0  sin£  costi 

(2.30b) 

03  *  cos0  sin<  sini) 

(2.30c) 

0*4=  sin0sin< 

(2.30d) 

03  -  sin0  cost) 

(2.30e) 

06  -  cos0  cost  sini) 

(2.300 

07  =  sinO  sini) 

(2  30g) 

08  -  COS0  COST) 

(2.30h) 

It  is  quite  evident  from  the  above  equations  htrv  the  coning  angle  0,  along  vrith  the 
orientation  angles  (T| ,t  ,♦  -  can  control  the  changes  in  each  orbital  parameter 

a.ej,  and  Q.  Judicious  choice  of  these  angles  can  lead  to  maximum  changes  in  orbital 
state  parameters . 


Special  Cases.  Jenkins  [Ref:  10)  derives  certain  conclusions  with  a  solar  sail 
system  which  is  in  one-to-one  resonance  and  vhose  motion  is  in  the  ecliptic  plane. 
That  analysis  shows  that  there  exist  unique  certain  sail  setting  providing  maximum 
changes  in  the  semi-major  axis  in  each  of  the  three  cases  below: 

a)  spinning 

b)  tumbling 

c)  coning 

By  singularly  varying  the  sail  setting  angles  in  the  perturbation  equations  (221  - 
2  30)  the  condition  for  maximum  change  in  semi-major  axis  can  be  graphically  shown. 
This  fact  is  later  shown  from  a  three  dimensional  perspective. 

Spinning  Case.  The  solar  sail  assumes  a  strictly  spinning  orientation 

when  the  coning  angle  9  equals  00  degrees.  There  is  no  "coning"  behavior.  This 
configuration  is  shown  in  Figure  2.11.  For  motion  in  the  ecliptic  plane  (i  -0)  with  the 

angular  momentum  vector  R  aligned  perpendicular  to  the  orbital  plane  (T|  -  C  -  0* ). 
this  “spinning"  orientation  affords  the  most  semi-major  axis  change.  As  will  be  shown 

later,  this  is  also  true  for  a  certain  range  of  orbital  inclinations  (0*  i  i  i  43* ).  At 
inclinations  greater  than  43* .  the  spinning  orientation  actually  results  in  a  decrease 
in  semi-major  axis. 


Fig.  2.11:  Spinning  Case 


I 


> 


I 


Tnmhiing  Cas*  This  is  the  case  where  the  angular  momentum  vector  H 

is  parallel  to  i  or  j  orbital  plane  reference  components.  With  H  parallel  to  the  i  orbital 
axis,  the  solar  sail  experiences  no  change  in  any  of  the  orbital  parameters.  This  is  the 
case  when  the  solar  sail  is  parallel  to  the  solar  radiation  and  hence  no  resultant  thrust 

component  is  produced.  Vith  fl  parallel  to  the  direction  vector,  the  changes  are  in  the 
radial  and  normal  direction.  This  results  in  changes  in  the  semi-major  axis  and  in  the 

inclination  angles  which  are  coupled  via  the  $  phase  angle.  The  strong 
implication  is  that  conditions  favoring  a  change  in  a  will  also  result  is  an 
accompanying  change  in  inclination  (which  can  be  unwanted). 

Coning  Case.  Coning  occurs,  when  the  sailcraft's  axis  of  symmetry  bj 

is  not  colinear  with  the  angular  momentum  vector  R.  Jenkins  investigated  this  coning 

motion  with  8  -  45*  and  discovered  the  importance  of  the  ($  -  x|>)  parameter  in 

controlling  the  magnitude  of  Aa.  and  that  there  exist  trade  off  possibility  between  Aa 

and  Ai.  the  resultant  change  is  inclination.  The  nature  of  this  coning  action  is  further 
investigated  in  this  study. 


% 


% 
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3-1  Dynamic  Sv«te«*  nprimi?*jinn 

SingJfcSlMfc  System.  Bryson  and  Ho  (Ref:  24)  describe  the  single-stage 
transition  from  some  initial  state  x(0)  to  a  near  state  x(l)  via  some  choice  of  control 
yector  u(0)  and  a  given  Operating  Function  F*  mathematically  as . 

x(  1 )  =  FMx(0),u(0)]  (3-D 

and  schematically  as .... 


Fig.  3-2:  Floar  Chart  for  a  Single-Stage  System 

(Ref:  24] 

vith  a  performance  index  of  the  form 


J  =  $  tx(  1 )]  =  L*  1 1(0),  u(0)]  (3  2) 

arhere 

x(0)  -  knoarn  state  at  initial  time  tg.  (n-dimensionai) 

u(0)  -  control  vector  (0 , 71 ,  £  ,  $  -  xjr)  (m-dimensional) 

x(  1 )  -  State  at  some  future  time  tj 
L*  (  1  -  Lagrangian  of  the  initial  state 


The  objective  of  an  optimization  problem  leads  to  the  maximization  of  this  performance 
index  J  .  Bryson  and  Ho  present  an  adjoined  performance  index  J  by  adding  (F*[l(0), 
ll(O)]  -  x(  1 )  -  0))  to  Equation  (3-2). 

3  .  ♦Ia(l)]*f  Ii(0),u(0)].Xt(1)(f1i(0),u(0)I-i(1)}  (3.3) 


vith  constant  Langrangian  multiplier  XT(  1 ) 

Introducing  H*  -  L*lx(0),  u(0)l  ♦  AT(l)F*[x(0),  u(0)]  (3.4) 

into  equation  (3-3)  results  in . 

J  -  ♦  Ix(  1 )  ♦  H*ll(0),  u(0),l(  1  )J  *  XTx(  1 )]  (3  3) 

The  total  derivative  of  J  is  expressed  as 

dj  =  (a»-  xT(n)  did)  +  air  dx(o)  „  air  du(o)  (36) 

dx(l)  dx(0)  du(0) 

With  the  judicious  choice  of  (1)  *  d$/dx(l).  Equation  (3  6)  simplifies  to  the 


foiloxring: 


w 


w 


where 

dH*/dx(0)  s  gradient  of  Jwrtx(O),  holding  u(0)  constant  and  satisfying  Equation  (31): 

and 

dH*/du(0)  a  gradient  of  J  wrt  u(0),  holding  x(0)  constant  and  satisfying  Equation  (31) 

Optimization  leads  to  finding  a  stationary  value  of  J  (i.e..  dJ  -  0)  for  a  given  initial  state 
x(0)  and  dx(0)  -  0.  Therefore,  the  stationary  value  of  J  is  found  if 

dJ  =  dH*/du(0)  =  0  (3.g) 


;•  Once  the  conditions  satisfying  Equation  (38)  are  met.  Jean  then  be  obtained.  Equation 

(3  8)  is  called  the  Optimality  Condition.  Applying  this  condition  to  Equation  (3.4)  this 
condition  states  that 


« 


4 


* 


* 


—  -  —  1 1(0),  u(0)  ]  .  —  [*V|I(0),U<0)]1  -  0 

du(0)  du  du 

Multipia-^tftga  System.  Progression  from  a  single-stage  system  to  a  multiple 
stage  system  can  be  easily  done  by  re-expressing  Equation  (3.1 )  as 

x(i  ♦  1 )  -  F* 1  x(i),  u(i)  ]  (3.10) 

where  the  initial  state  x(0)  is  given  and  i  -  0,  1.  ...  N-l  with  N  s  nth  stage.  Here. 
Equation  (3. 10)  describes  a  sequential  set  of  equality  constraints  with  x(i)  as  a 
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sequence  of  n-vectors  to  be  determined  by  another  sequence  u(i)  of  m-vectors 
Schematically,  this  is  expressed  as  a  cascade  system. 


(Ref:  24] 

The  corresponding  adjoined  performance  index  J  is 

3  -  ♦  li(n)l  L‘ [i(i),  u(i)]+  AT(i*l){FMA(i),u(i)  -A(i*l)}].  (3.11) 

Using  similar  operation  as  for  the  single-stage  case,  the  optimality  condition  is 
expressed  as 

d^/duti)-  0  (312) 

vhere  H*  =  L*  I  l(i),u(i)  ]  +  AT(  i  *1)  {fM  l(i),u(i)  ]  -  l(i+l)}  (313) 

for  i  -0,1, ...  N-l 

Hence,  to  find  a  control  vector  sequence  u(i)  that  produces  a  stationary  value  of  the 

performance  index  J.  effort  must  be  placed  in  solving  the  following  difference 
equations: 

x(i  +  1)  =  F1  (  x(i),u(i)  ]  (3  14a) 

and 


X(i)  =  IdFVdxd)]1  A(i  ♦  1 )  ♦  IdLVdxd)]1 


(314b) 


u(l )  is  that  control  sequence  that  males  flotation  ary  (From  Equation  3  12)  That  is 


dHV  du(i)  =  dL‘  /du(i)  ♦  AT(i  ♦  1 )  dF*  /du(i)  *  0  (315) 

(REF:  24] 

Application  to  Problem.  The  single-stage  optimization  can  easily  be  applied  to 
the  present  objective  of  maximizing  the  final  semi-major  axis.  The  problem  parameters 
▼ere  established  as  follows: 


let  U(0)  -  (0,7), <,♦-♦)  -  (8.7), CPA)  (316a) 

x(0)  =  initial  state  (3 16b) 

i(l)  =  final  state  (316c) 

i  =  [  a ,  e  ,  i  ,Q] 

F*  [  1  =  lxo  +  Ax]  (3  16d) 

A  =  I  1  (316e) 

L*  =  0  (3 16f) 


Applying  equation  (3-16)  to  the  single-stage  system  optimality  condition,  the  condition 
to  be  met  is  found  to  be 


dH7du(0) 


AT(1)  dFa/du(l)  =  0 


(317) 


For  maximiziflg  the  semi-major  axis,  the  Lea  grange  multiplier  vector  becomes 


A11'  -  11,0,0,01 


(3.18) 


Since  F*  constitutes  the  changes  in  the  state  vector,  we  can  express  it  formally  as.. 


F*  - 


V  Aa 

e#  ♦  Ae 


i.  ♦  Ai 
Q„*  AQ 


(319) 


Similarly,  the  control  vector  u(0)  is  expressed  as 


lOI,01.C,,<8**)ll 

>ul'.  u22.  “31.  u<‘1 
(Uj1)  for  i  -  12.3.4 


(320) 


By  definition,  the  expression  dF*/3u  ^  becomes  (with  PA  -  $+\J/), 


dF’/dU1 


dFg/dfl  dFa/dTl  3Fa/d<  dFa/3PA 
aFe/ae  aFe/aTi  aFe/a<  dFe/aPA 


aFj/ae  aFj/a-n  aF/a<  aF/aPA 
aF Q/a0  aFQ/aTi  aFQ/ac  aFQ/aPA 


( 3.2 1 A ) 
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Since  XT(1)  is  not  identically  zero  in  equation  (317),  the  following  must  hold  true. 

dFa/du1  =  1C]  =  10].  (3.22A) 

The  scope  of  the  problem  is  trmendously  reduced  by  setting  the  performance  index  to 

be  the  optimum  change  is  Aa  per  orbit.  Interest  then  lies  only  in  the  function  Fa 
which  constitues  the  initial  semi-major  axis  and  its  change.  Furthermore,  equation 
(3  21  reduces  the  complication  one  more  degree  by  taking  the  partial  derivatives.  Since 
a0  is  the  initial  state  and  is  constant,  its  presence  in  the  differentiation  can  be  ignored. 

The  remainder  is  plain  Aa.  The  scope  of  the  problem  in  hence  greatly  reduced  by  an 
order  of  magnitude 

This  means  that  equation  (3  21A)  can  be  reduced  in  size  to  the  following. 

dFa/du(l)  -  [  dFa/d8,  dFa/dn,  dFa/d<,  dFa/dPA]  (3  21B) 

using  the  notation  expressed  in  equation  (3  20)  this  can  be  rewritten  as 

aFa/dud)  - 1  aFa/dUj,  aFa/au2,  aFa/au3  aFa/a  u3i  (3.210 

or 

aFa/au(o  -  [aFa/auji  for  i  - 1.2.3. <  (3.21D) 

The  control  vector  u(l)  that  satisfied  this  Optimality  Condition  Equation  (322)  will 
determine  the  stationary  value  of  the  performance  index  J.  where  J  is  stated  above  as 
the  maximum  final  orbit.  J  -  a(N).  For  this  single-stage.  J  -  ad)  -  the  semi-major 
axis  after  1  orbit.  The  search  for  this  control  vector  is  the  quest  of  this  effort.  The 
solution  approach  follows. 


(• 

I* 


,t 


,  4 
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3  2  Solution  Approach. 

Solution  Formulation.  F*  is  given  as  a  function  of  the  control  parameter  (up. 
Specifically,  looking  at  Fa  only. 

Fa  =  Fa(u,,u2,  u3,  l»4)  =  FJup  for  i  -  1. 2.3.4  (3.23) 

where  (u,,  U2.  U3,  U4)  -  (0 , 71 .  C  ,  8"*)  (3  24) 

Fft  is  essentially  a  single  function  of  four  variables.  The  interest  is  in  finding  values 
of  (Uj)  at  which  Fa  is  a  maximum.  Under  certain  conditions,  at  such  "places"  in  the 
space  of  variables  (up. 

[C]  -  [aFj/dUj]  -  10]  for  i  -  1 ,2,3.4  (3  23) 

One  method  of  searching  for  this  special  (up  is  as  follows: 

a.  assume  an  initial  guess  u’  and  that  [  dFa/dllj*]  exist. 

b.  assume  that  Fa(  up  -0  can  be  expanded  as 

F^Uj.U^Uj.U^)  =  Fa(ul,,U2*,U3*,U4*) 

+  1  3Fa/dUj  (up.Uj'.Uj'.u^*)  (Uj-  Uj’) 

*  *  I.  l,JFa/(aujau.)  (u/.Uj’.Uj-.u,-)  (ur  Uj')  (u„-  uB')...  (3  26) 

and  dFa/dui(u,,u2,u3,u4)  =  dF^a^  (uI*>u2*,u3*,u4*) 

I  d2Fa/(duiduj)  d2Fa/(dujaum)  (up.u^.Uj'.u^)  (uj-  up 


(3  27) 


[  PCjj]  =  ^/Oujdu,)  (Uk‘)  =  dC/OujU.)  (Uk*) 


(328) 


in  Equation  (3-27)  and  using  Equation  (3.2?),  gives 


Vui‘  -  |  IKijl-'ICl  (ul,.u2‘.u3-.u4*). 


(3.29) 


Defining  5u  =  (u,-  l», *). 


(3.30) 


Equation  (329)  can  be  expressed  as 


Su  -  -£  |PC]  ij-'tci  (u .•). 

>■1  1 


(3.31) 


Understanding  that  the  summation  must  be  taken,  a  simplified  expression  can  be  stated 


Su  -  -lac/aui’1  IC1  (up. 


(3  32) 


In  Equation  (3  32).  [C]  (u.* )  denotes  the  Optimality  Vector  evaluated  at  the  initial  guess 


control  vector  U.*  fori  -  1. 2.3.4. 


Su  denotes  the  required  change  in  control  vector  u  that  would  satisfy  the  condition 
for  a  stationary  performance  index  J.  Iterations  of  an  "initial  guess"  is  necessary  until 

this  Su  ~  0  .  A  reduction  of  the  iteration  frequency  can  be  achieved  by  providing  a 
“good"  initial  guess.  A  graphical  method  of  determining  this  guess  is  presented  in 
Section  IV  under  Surface  Representation. 
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Analytical  VS  Numerical.  Solution  to  Equation  (3  32)  requires  the  evaluation  of 
the  partial  derivatives  of  four  state  parameter  functions  (Fa,  Fe.Fi,  F#)  with  respect 

to  four  control  variables  ($,  T|  ,£  ,PA).  This  involves  the  determination  of  a  [4x4] 
matrix  and  its  derivative.  Recall,  the  objective  is  to  find  the  control  parameters 
leading  to  a  maximum  change  in  the  semi-major.  From  the  C  Equation  (3  21).  it  is 
apparent  that  interest  lies  in  the  1st  row  elements.  This  row  constitutes  the  optimality 
vector.  This  reduces  the  main  effort  considerably.  This  simplification  provides  the 
following: 

[C]  =  [dAa/d  6  dAa/dT]  dAa/d<  dAa/dPA]  (3  33) 

Disguised  in  Aa  are  horrendous  amounts  of  differentiation  that  can  easily  lead  to  any 
amount  of  errors.  The  choice  made  here  was  to  circumvent  the  series  of 
differentiations  and  use  numerical  techniques  to  determine  the  elements  in  Equation 
(3  26).  To  appreciate  this  approach,  one  must  initiate  the  series  of  differentiation. 
Appendix  B  provides  the  reader  this  insight.  The  reason  for  the  numerical  method 

choice  and  why  the  analytical  approach  is  avoided  is  apparent. 


3  3  Numer ical  Fgrmyll;Hi9P 


The  Numerical  Models.  As  a  first  approximation  to  the  partial  derivative  of  the 
multi-variable  function  F(ui),  the  numerical  differentiation  techniques  called 
Newton's  Forward  and  Backward  were  used.  These  are  expressed  as 

F(U|  +  Suj)  -  F(U|) 

F'tUj)  =  _  (Forward)  (3.34) 

SUj 

F(Uj)  -  F(Uj  +  SUj) 

and  F'(Uj)  *  _  (Backward)  (335) 


For  a  simple  two-variable  case,  these  can  be  interpreted  as  the  slopes  of  the  two  lines 
shown  in  Figure  3  8. 


Fig.  38  Derivative  Approximation  Models 


Initial  Results.  Initial  Results  were  not  encouraging.  Initial  intuitive 
expectation  vere  filled  with  optimism  for  a  semi-veil  behaved  derivative  of  the 
function  F(Uj).  This  was  not  the  case.  These  formalas  displayed  a  "wandering"  type 

behavior  with  the  approximated  derivatives  for  each  iteration.  The  Forward  Formula 
can  be  re-expresssed  as 

Au  «  {F(u)new  -  F(u)old)/8u  «  -  F(u)/F’(u)  (3  46) 

Using  this  formulation,  the  expectation  was  for  some  kind  of  rough  convergence. 
Instead,  the  Au^  were  fluctuating  after  each  updated  estimate  F(Uj)nev. 

Sources  of  Errors.  A  study  of  test  cases  suggested  that  approximate 
derivatives  obtained  from  from  such  polynomial  F(Uj)  be  viewed  with  skepticism  unless 

very  accurate  data  are  available.  Even  with  accurate  data  (or  initial  guess  as  it  is 
commonly  called),  the  accuracy  diminishes  with  increasing  order  of  the  derivatives. 
This  is  the  problem  dealt  with  the  second  order  derivatives  .  The  dominant  error 
source  is  in  the  input  errors.  These  proved  very  critical.  Even  when  the  initial  guess 
was  close  to  the  theoretical  value,  the  input  was  still  very  critical  because  the 
approximating  (along  with  the  perturbation)  algorithms  magnify  them  enormously. 

The  crucial  factor  seemed  to  be  the  magnitude  of  Su.  The  magnification  of  input  error 
behaved  inversely  to  this  value  whereas  the  inevitable  truncation  error  was  directly 
affected. 

Improved  Model.  The  models  intially  used  were  abandoned  after 
countless  futile  attempts  to  contain  the  "wandering”  derivatives.  From  Figure  3-8  it  was 
evident  that  a  more  accurate  approximation  would  be  achieved  by  using  the  slope  of  a 
line  connecting  points  A  and  B.  This  approach  is  commonly  called  the  Stirling  Method. 


It  is  basically  the  Newton  method  modified.  This  veil-known  formula  also  known  as 
central  differencing  is  expressed  as 

F'(u)  *  {F(u  +  8u)  -  F(u  -  Su)}  /(28u)  (Stirling)  (3.47) 


Numerical  Formulation.  This  formula  was  employed  in  determining  the 
required  derivative  expressed  in  Equations  (3  44)  of  the  previous  section.  Appendix  B 
shows  how  it  was  coupled  to  accomodate  a  function  of  four  variables. 

Results.  The  use  of  Stirlings  Method  provided  better  results.  The 
"wandering"  phenomenon  experienced  with  the  previous  Newton  s  methods  was 
greatly  diminished  by  an  order  of  magnitude. 

From  the  perspectives  peak  points  were  selected  as  test  cases.  Known  maxima  is 
inputed.  With  the  appropriate  coordinates  for  the  maxima  (from  perspectives  as  input, 
convergence  to  within.  .1  degree  was  achieved.  This  related  the  fact  that  a  maximum 
did  indeed  exist  within  the  area  of  search.  This  gave  additional  credence  to  the 
perspectives.  The  various  identifiable  coordinates  were  checked  with  similar  results. 
Convergence  was  achieved  in  two  iterations.  The  bottom  line  was  that  the  required 
calculated  change  in  control  parameters  was  within  .1  degree  a  change  so  small  that 
virtually  a  maximum  was  discovered.  This  was  satisfaction  but  not  with  surprise  since 
approximate  oordinates  of  the  maxima  were  inputted.  The  real  test  follows. 
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Annrntimiiftg  32F/^ui^ui)  The  Stir  ling  Method  is  used  to  approximate 
the  d2F/(dujduj)  matrix.  For  simplicity,  the  following  definition  is  made: 

Let  (PCi  j)  =  &Y/  dUjdUi  -  ac/  aUj  (3.36A) 

where  C  -  dF/  dUj  (3  36B) 

The  subscripts  i  and  j  indicate  a  particular  control  vector  from  the  set  (u  j.  u2,  u^,  u^)  - 
(0.  71,  £,  PA).  Recall,  that  for  an  initial  guess  control  vector,  F  is  defined  as 

F  =  Ffllj*,  U2*.  U3*.  U4‘)  (3.37) 

where  the  superscript  (*)  denotes  the  initial  guess.  Note  that  for  if  any  iteration  is 
desired,  the  Uj*  term  would  be  the  updated  value  which  would  become  the  new  initial 
guess. 

Hence.  (Uj\  u2\  u^’,  u^’)  -  (0*.  T|“ ,  £*.  PA’)  *  initial  guess. 

Vith  this  notation,  equation  (3  36)  can  easily  be  expanded  using  Stirlings  Method. 

pFamnift  To  show  how  the  matrix  is  constructed,  the  PC)2  component  is 
selected  as  an  example. 

>2h 

PCU  -  -  (3  38) 

dujdu2 
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By  following  the  same  procedure  for  the  other  elements,  I  PCjj]  matrix  can  be 
constructed  as  follows: 


d2Fa/(du,)2 

d2Fa/(du2duj) 

32Fa/OujdU|) 

d2Fa/(du4du,) 


d2Fa/(du,du2) 

d2Fa/(du2)2 

d2Fa/Oujdu2) 

d2Fa/(du4du2) 


d2Fa/(du,du5) 

d2Fa/(du2du3) 

d2Fa/(du3)2 

d2Fa/(du<du3) 


d2Fa/(du,du,|) 

d2Fa/(du2du4) 

d2Fa/(du3du4) 

d2Fa/(du4)2 


(3.44) 

This  matrix  is  used  in  computing  the  change  in  control  parameters  necessary  to 
determine  the  stationary  value  of  the  performance  index  J.  From  Equation  (3  ). 

Aui  *  “new  -  “old  "  '  <c*  V  IC1  for  1  * 12.3.4  (3-43A) 


where  (C* }  denotes  the  dFa/duj  evaluated  at  the  initial  control  values  and  [C]  denotes 
the  d2Fa/(dUjduj)  -  (PCjj].  The  C  notation  is  later  used  in  the  algorithm  developed  to 
compute  the  Auj.  Hence, 


A“i  '  “new  -  “old  '  -  <c*  V  tPCijl  (3  45B) 

Note  that  (C* )  is  a  four-element  row  vector  and  [PCij  ]  is  a  (4x4)  matrix.  Performing  the 
appropriate  operation,  the  required  change  in  the  control  vector  sought  here  is 
plainly. 


A“i  ’  “new 


-“old 


-  ‘PCij)'1 


(C*) 
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(3.450 


A  stationary  value  is  achievable  if  there  exist  a  control  variable  set  Uj  such  that  the 
Auj  is  approximately  zero  That  is. 

Aui  “  unew  '“old  -  ‘  (pcijrI  “  (0)  (3  43D) 

This  vould  indicate  that  there  are  no  changes  necessary  in  the  previous  control 

variables  set  to  arrive  at  the  optimum  Aa.  Such  is  the  particular  1  set  of  controls  that 
vould  render  the  performance  index  J  stationary. 

Software  Development.  Equation  (3  45D)  is  the  basis  of  the  software  that  is 
developed  to  find  the  stationary  value  of  the  Performance  Index  J.  The  process 
involves  the  coding  of  the  perturbation  solution  equations  presented  by  Jenkins 

[Ref:10l.  The  Aa  equation  is  the  only  equation  generating  the  required  vectors  and 
matrices  in  Equation  (3-4)  since  the  objective  solely  called  for  maximizing  this 
function.  The  incorporation  of  the  other  state  functions,  as  expressed  in  Equation 
(3.19).  can  readily  be  employed  with  little  transitional  difficulty.  The  mechanism  of 
the  search  process  included  several  iterations  to  construct  the  elements  of  (C)  vector 
and  [PC]  matrix.  Once  this  [PC]  matrix  was  found,  a  specially  formatted  program  is  used 
to  compute  its  inverse  PCI.  Equation  (3.45D)  is  then  incorporated  to  determine  the 
neccesary  changes  in  the  control  vector  u.  The  process  continues  until  a  specified 
convergence  criterion  is  established.  The  interactive  program  outputs  all  the 

computed  (intermediate  and  final)  vectors  and  matrices  along  with  the  required  AU; 


3  4  The  Nature  Of  The  At  Function. 

The  change  in  the  semi-major  axis  is  the  main  objective  of  this  study.  It  almost  seems 
as  a  trivial  task  until  the  function  itself  is  confronted.  Aa  is  a  function  of  four 

variables.  Aa  =  Aa(0,T),C  The  main  difficulty  in  maximizing  this 

function  analytically  is  with  the  complexity  of  the  derivative  of  the  function  (as  was 
explained  in  section  3-1  and  Appendix  B).  The  vise  choice  of  employing  numerical 

techniques  to  determine  the  maxima  of  Aa  still  presents  another  difficulty-the  initial 
guess.  This  initial  quess  becomes  the  starting  point  in  the  search  of  the  set  of  control 

vector  which  maximizes  Aa.  Due  to  the  search  pattern  of  the  modified 
Newton-Rhapson  (or  Stirling  Method),  this  initial  guess  becomes  critical  and  must 
serve  as  a  "good"  starting  point.  This  is  where  the  difficulty  begins:  What  is  a  "good 
initial  guess"?  How  can  a  region  of  "good  guesses  be  established?"  The  answers  can 

be  provided  by  a  surface  representation  of  the  function  Aa.  The  next  section 
explains  an  In-House  capability  available  through  UNIX/VAX  system  that  provides 
this  surface  representation. 


3.5  Surface  Representation 


A  canned  plotting  package  called  "S-Package"  residing  in  the  UN IX/ VAX  system  is 

used  to  provide  three-dimensional  perspectives  of  the  Aa  function.  This  program 
provides  an  quasi-isometric  three-dimensional  perspective  of  an  output  function 
with  respect  to  two  input  variables  with  the  output  function  as  the  third  dimension. 

Analyzing  Aa  is  accomplished  by  holding  the  angular  momentum  orientation  angles 

(T|  and  O  stationary  while  allowing  variations  in  the  coning  angle  (0)  and  the  phase 

angle  (PA  -  The  resulting  perspectives  show  the  behavior  of  the  Aa  function 
as  these  selected  parameters  are  varied.  These  are  presented  in  more  detail  in  the 
Results  section.  Such  surface  representations  describing  the  behavior  of  the  change 
in  semi- major  axis  with  respect  to  a  coning  Solar  Sail  have  never  been  seen  before! 


These  perspectives,  provide  very  valuable  information  on  the  nature  of  the  Aa 
function  as  certain  orientations  are  approached.  More  importantly  the 
perspectives  serve  as  a  source  of  information  on  the  region  of  search  from  which 
one  can  determine  a  "good  initial  guess".  The  maxima  and  mimina  are  apparent. 
Additionally,  problem  areas  that  might  cause  a  gradient  search  scheme  to  zander 
excessively  and  perhaps  fail  can  be  avoided  simply  by  identifying  critical  areas  from 
the  perspectives. 

No  matter  what  scheme  is  used  to  present  a  function  of  more  than  two.  there  still 
remains  the  question  of  the  behavior  of  all  the  variables.  Man.  limited  to  working 
with  three  dimensions,  representing  a  function  of  more  than  three  variables  is  quite 

a  task.  In  this  study,  Aa  is  a  function  of  four  variables;  therefore,  a  five-dimensional 
perspective  is  ideal  to  explicitly  show  its  true  behavior  with  respect  to  all  the 
variables.  Accepting  this  limitation  and  expressing  the  function  with  respect  to 
three-dimensions  is  the  best  anyone  can  presently  do.  Despite  the  limitation,  the 
benefits  of  having  at  least  a  five-dimensional  representation  is  overwhelming. 


IV.  Results 


41  Surface  Perspectives 

The  nature  of  the  Aa  function  is  analyzed  for  the  occurrence  of  maxima  as  the  various 

control  parameters  (8,  T).  £,  PA)  are  varied.  Being  a  multi-variable  function,  it 

becomes  practically  impossible  to  loot  at  the  overall  nature  of  this  Aa  function  as  all 
the  parameters  are  allowed  to  vary.  It  is  instructive  to  hold  any  two  less  dynamic 
variables  and  perspective^  look  at  the  function  with  the  two  more  dynamic  variables 

allowed  to  go  free.  This  would  constitute  a  three-dimensional  representation  of  the  Aa 
function.  Abundant  information  can  be  read  off  such  perspectives  since  all  the 
variations  and  and  resultant  behavior  of  the  function  can  be  readily  seen  in  a  given 
perspective.  Anomalies  are  conspicous  by  their  presence  (if  any  are  present).  The 
extrema  (maxima  and  minima)  can  easily  be  spotted  and  located  for  the  conditions 
given.  It  is  also  very  helpful  in  identifying  areas  of  "stagnation"  or  "stability"  in 
which  the  function  can  very  well  remain  stationary  even  with  small  perturbations. 
This  very  fact  that  one  has  a  three-dimensional  interactive  look  at  the  behavior  of  the 
function  makes  a  valuable  preliminary  tool  in  the  optimization  process  where 
attaining  a  "good  initial  guess"  is  so  paramount.  Such  information  is  available  with 
perspectives.  Such  dynamic  information  is  found  absent  in  tabular  formatted  data. 

The  information  contained  in  the  perspectives  found  in  this  study  contributes  to  the 
search  scheme  developed  to  find  the  optimal  control  setting  angles  to  achieve  an 
optimum  change  in  the  semi-major  axis  per  orbit.  The  "initial  guess"  barrier  has  just 
been  broken. 

Variation  &  Inclination.  Interest  exists  in  the  behavior  of  the  change  in  the 
semi-major  axis  Aa  due  to  changes  in  the  sailcraft's  orbit  inclination.  Figure  4.1  shows 
the  relative  Aa  magnitude  of  three  select  inclination  with  changing  coning  angles. 
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Note  the  decreasing  Aa  magnitude  as  the  inclination  increases  and  then  its  subsequent 
increase  to  form  two  maxima  as  the  inclination  is  greater  than  45* .  This  behavior  was 
graphically  evaluated  more  closely  as  inclination  varied  from  0  to  90  degrees.  Figures 
4.2 A  through  42G  (presented  in  the  following  pages)  show  the  peculiar  behavior  of 

this  Aa  function  as  the  inclination  is  increased.  Each  figure  depicts  a  different  case 

vith  with  its  corresponding  Aa  magnitude.  Although  these  relative  Aa  magnitude 

differences  are  not  shown  in  these  figures,  the  maximum  Aa  experienced  per  case  is 
indicated.  The  important  information  is  in  the  overall  behavior  of  this  change  in 

semi-major  axis.  The  interesting  aspects  of  this  Aa  behavior  are  summed  as  follows: 


As  the  inclination  is  increased  from  0*  to  90* , 
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a.  the  magnitude  of  Aa  decreases  as  the  inclination  approaches  45 
degrees. 

b.  the  magnitude  of  Aa  increases  as  the  inclination  becomes  greater 

than  43  degrees  but  does  not  attain  its  original  maximum  value. 

c.  the  local  maximum  shifts  location  away  from  0  -  90*  (at  i  -  0’ )  to 
form  two  maxima  at  8j-35‘ and  02- 143 ‘  (ati-90’); 

d.  maxima  occur  at  PA  --90‘  and  remains  stationary  at  this  value,  i.e., 
it  does  notvary  with  changes  in  inclination; 

The  three  select  cases  are  further  examined  in  the  following  pages  u>  show  why  such 
particular  control  angles  work  to  maximize  Aa.  The  three  cases  looked  at  are  as  follows; 

Case  1;  0*  inclination 

Case  2:  43'  inclination 


Case  3;  90*  inclination 


id4  DU 


Max  Ao  =  1.7530  x1040U 
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Max  Aa 


Max  Aa  =1 .70097  x  10  DU 


Max  Aa  =1 .70097  x  10  DU 


4 


4 


% 
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0  Degree  Inclination  Motion  in  the  ecliptic  plane  with  the  coning 
angle  set  at  90*  and  the  phase  angle  at  -90*  has  produced  the  largest  change  in  the 
semi- major  axis  at  2.337 xlO'*  DU.  Figures  4.1  compares  the  behavior  in  this  inclination 

with  other  inclined  orbits;  Figure  4.2A  perspective^  portrays  the  behavior  of  the  Aa 

function  as  any  of  the  other  control  parameters  (B  and  PA)  are  varied.  For  a  coning 
angle  of  90*.  the  solar  sail  is  (at  certain  intervals  in  the  orbit)  perpendicular  to  the 
solar  radiation;  in  such  an  orientation,  there  is  maximum  thrust  produced  in  the 
direction  away  from  the  sun.  However,  maximum  thrust  does  not  always  equate  to 

maximum  Aa.  It  could  very  well  work  against  maximizing  Aa;  such  a  case  is  when  the 
thrust  vector  has  a  component  working  against  the  orbital  velocity  vector.  Correct 
phasing  angle  PA  allows  the  optimum  use  of  maximum  thrust  conditions  leading 

towards  maximizing  Aa. 


This  PA  parameter  relates  the  initial  point  of  the  sailcraft  in  the  orbit  and  its  initial 
orientation  wrt  the  Body-centric  Inertial  IJK  Frame.  For  0  -  90* ,  the  latter  information 
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(proper  phase  angle)  is  inconsequential  since  the  sail  is  in  a  pure  spin  motion. 
However,  for  any  other  non-pure  spin  alttitude.  this  information  is  of  critical 
importance.  Hence,  angle  PA  only  relates  where  the  sailcraft  is  in  the  orbit  for  this 
pure  spin  configuration.  At  PA-  -90*.  the  sailcraft  is  configured  as  shown  above  in 
Figure  4.3;  Note  where  the  sailcraft  is  allowed  to  start  its  journey  from;  90* 

counter-clockwise  from  the  i  axis.  This  orbit  configuration  produces  the  most  Aa(Aa 
-  2.337  xl(T4Du)  the  torque-free  solar  sail  can  muster  optimally  at  improving  its  radial 
distance  from  the  planet  earth. 

H  Degree  Inclination.  As  seen  earlier.  Figure  4.1  provided  an 
indication  of  what  would  happen  to  Aa  when  the  orbit  inclination  is  increased.  The 
shift  of  the  single  maximum  at  (8,  PA)  -  (90,-90)  to  form  a  pair  of  maxima  can  be  very 

clearly  seen  in  the  perspective.  Figure  4.2D.  No  longer  does  0  at  90  degrees  monopolize 
the  maximum  point  at  this  inclination:  the  two  maxima  occur  at  (33*.  -90* )  and  (123*. 
-90* ).  This  is  important  to  understand:  the  pure  spin  configuration,  though  it  may 
present  the  most  exposed  surface,  does  not  present  the  optimum  condition. 


Fig.  4.4A:  43*  Inclination 


Figure  4.4A  above  shows  the  various  positions  in  the  orbit  relation  to  the  attitude  of  the 
solar  sail  wrt  the  sun.  Only  the  35*  coning  angle  case  is  depicted  here.  The  125* 
coning  case  contributes  the  same  effect.  This  is  studied  next. 

The  relationship  between  the  two  coning  options  are  investigated  to  determine  the 
resulting  components  of  thrust  along  the  velocity  vector.  Both  diagrams  shown  here 
in  Figure  4.4B  relate  to  the  position  of  the  sailcraft  at  its  starting  position.  Since  the 
sailcraft's  precession  rate  is  equated  to  its  orbital  mean  motion,  the  coning  motion  will 
be  cyclic  and  in  phase  with  the  orbital  period. 


From  Diagram  A: 


From  Diagram  B 


Ty  -  T  cos  33*  -  819  T 

Tv  -  T  sin  33*  -  819 T 


Tv(0  -33*)-Ty(0  -123*) 


Since  the  two  components  are  identically  equal,  their  corresponding  da  effects  will 
also  be  equal.  For  the  configurations  of  higher  inclinations,  this  offers  an  option  in 
the  choice  of  coning  angle. 


9&  Degree  Inclination-  This  specie!  inclination  is  depicted  in  Figure 
(4.3)  below  Fimpie  (Ref:  31  states  that  this  orientation  affords  the  sailcraft  continuous 

energy  increase  throughout  the  orbit.  However,  the  magnitude  of  Aa  is  smaller  than 
for  that  seen  in  the  0*  inclination  case  because  of  the  small  thrust  component  in  the 
direction  of  the  orbital  velocity  vector. 


velocity  vector  (see  Section  II.  Previous  Efforts).  Given  the  similar  orientation  Garvin 

refers  to,  this  tilt  angle  is  just  the  coning  angle  8  of  this  study.  Tsu  (Ref:  19]  also 
arrives  at  this  optimum  till  angle  of  33  degrees.  Configuration  B  depicts  this  control 
setting.  For  comparison  vith  the  same  setting  but  at  a  different  phase  angle  (PA  - 

-270* ).  the  thrust  component  is  no  longer  vorking  to  Aa.  but  vorking  to 

minimize  Aa.  In  fact,  the  minimum  Aa  is  experienced  vith  this  control  setting. 


Variation  in  PA.  If  the  variation  in  8  is  restricted  0  <  8  <  360  *.  a  maximum  Aa 
occurs  at  PA  equal  -90*.  Such  occurence  at  this  value  can  be  seen  analytically  by 
studying  Equations  (2.21)  and  (2.22)  shown  below: 


Aa  =  a  L5  DTP/m  [da,  +  da2  +  da3  +  da^  +  da5+  da6]  (2.21) 

da,  =  '/4D,  2[(p2+3p4)cos(^)-(3pl+p5)9in(^-t)l  (2  22a) 

da2  =  Vx  D,  D2  I(p2-(H)  8in($-\jr)  +  (p5-pl  )  cos($-*)  1  (Z2Zb) 

da3  *  2  D,  D3  (  06  cos($->|/)  +  03  sin($-^)  1  (2.22c) 

da4  =  2  D2  D3  (  £3  cos(4>-^)  +  p6  sin($-*) )  (2.22d) 

da3  =  14D2D2  l  (3p2+p4)  cos(<|>-*)  -  (pi  +3p5)9m(*-4f)  I  (2.22e) 

da6  *  D32  l  (02 +p4)cos($-^)  -  (01  ♦  05)  sin($-*)  1  (2.22f) 


From  Equation  (2.22).  one  can  say  that  maximizing  the  summation  of  the  da  terms 
equates  to  maximizing  Aa.  A  plot  (see  Figure  4.6)  of  these  da  terms  above  for  variations 

in  the  PA  (recall.  PA  -  $-xJ/)  component  explicitly  shows  that  the  summation  is  indeed 
maximum  at  PA  equal  to  -90  degrees. 


Variation  ifcg  Angular  Momentum  Vector  flciglUatiQfl  The  angles  Tl  and  £  which 

describe  the  angular  orientation  of  the  angular  momentum  vector  fl  with  respect  to 
the  orbital  reference  frame  (ijk)  has  been  superficially  neglected.  In  this  study,  they 
have  taken  on  common  values  of  0*  or  90* .  Earlier,  these  angles  were  termed  the 
“lesser  dynamic'"'  parameters  and,  therefore,  vere  not  varied  in  the  previous 

perspectives.  On  the  contrary,  variations  is  Tl  and  C  have  sufficient  control  on  the 

behavior  of  Aa  or  else  they  would  not  have  been  included.  Their  major  is  the  shifting 
the  location  of  the  maxima.  Here  are  sample  variations: 

Fixed  Tl.  Figures  4.7A.B.C  show  just  such  behavior  as  T|  was  held  fixed  at 
zero  and  4  was  allow  to  vary  from  30*  to  90*  for  motion  in  the  30*  inclination  plane.  As 

£  is  increased,  the  maxima  no  longer  becomes  stationary  at  PA  -  -90  but  seems  to 
propagate  or  shift  towards  PA  -  33*  ■ 

Fixed  C.  A  similar  but  out  of  phase  situation  occurs  when  (  is  fixed  at 

zero  with  Tl  allowed  to  vary.  At  Tl  -  30* .  the  Aa  function  assumes  the  same  shape  as  in 

the(Tl,  O  -  (60.0).  This  case  is  depicted  by  Figures  4.8A.B  and  C.  The(Tl.C)  ■  (90,0)  is 
essentially  the  tumbling  case  Jenkins  (Ref:  10]  investigated  earlier. 


Max  Aa  =  1.606 


I  =  ovkou 


Perspective  (i  =30*) 


2  1912 


$£*£  I  =  ovx©u 
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The  search  of  the  optimum  set  of  controls  to  achieve  a  maximu  change  in  a  per  orbit 
required  the  use  of  the  perspectives  as  a  guide  in  establishing  a  good  intial  guess  to 
start  up  the  search  algorithm.  After  noting  the  location  of  the  maiima  from  the 
particular  perspective  of  interest,  convergence  is  almost  expected  at  that  graphical 
coordinates. 

Observations.  In  employing  the  Stirling  Method  in  search  of  the  required  set  of 
optimal  control  vector,  a  few  peculiar  observations  were  made..  The  Stirling  Method 

uses  a  step  size  Su  to  increment  through  the  local  area  of  interest  in  search  for  the 
zero  of  the  function. 

a.  The  search  scheme  is  very  sensitive  to  the  step  size  Su.  Vhen  the 
coordinates  of  the  MMim*  (from  perspectives)  are  inputted,  the  value  of  the  (C* ) 
vector  components  become  very  small  often  in  the  neighborhood  of  1  x  10  "7  .  Since 

(C*)  is  the  initial  dFft/du;  evaluted  at  the  initial  input  control  values,  using 

near-maxima  control  values  would  make  the  incremental  step  size  too  small  to 
differentiate  a  change  that  can  be  discernable  by  the  search  scheme.  This  is  where 
larger  step  sizes  can  be  used  without  losing  any  appreciable  accuracy. 

b.  If  the  conditions  generating  the  perspectives  are  used  as  initial  guess  inputs, 
the  maxima  depicted  can  be  reaccomplished  to  a  higher  degree  of  accuracy.  Test  cases 
show  this  accuracy  to  be  good  to  the  seventh  order. 

c.  The  degree  of  accuracy  beyond  the  seventh  order  seems  to  be  hampered  by 
round-off  or  truncation  error  despite  the  fact  that  all  computations  are  accomplished 
in  double-precision  mode. 


Test  Cases.  Several  test  cases  were  made  primarily  with  the  use  of  the  data 
extracted  from  Figures  2.4A  through  2.4G.  Since  the  intent  of  the  preliminary  searches 


is  to  use  these  figures  as  the  primary  feedback  for  the  behavior  of  the  8a  function  and 
the  verification  of  the  output  of  the  search  algorithm,  deviations  from  the  initial 
conditions  generating  the  figures  vere  not  made.  The  results  on  these  preliminary 
evaluations  can  be  summed  up  as  follows: 

1.  Searches  using  coordinates  of  maxima  from  the  perspectives  yielded  similar 

Su  required  essentially  equal  to  (0.0.  0.0.  0.0,  0.0).  This  implied  that  the  point  of 
interest  is  a  maxima  to  begin  with  and  that  no  changes  in  the  control  parameters  are 
required.  Convergence  is  achieved  in  a  single  iteration.  Further  forced  iterations 
resulted  in  the  same;  i.e.  when  search  was  continued  until  an  iteration  limit  was 

reached,  the  required  Su  remained  at  zero.  This  displayed  stability.  Larger  step  sizes 
did  not  affect  the  convergence. 

2.  Deviations  in  any  single  control  parameter  (holding  the  other  three  at  initial 
maxima  condition  values)  resulted  in  more  than  one  but  less  than  three  iterations  to 
converge  to  a  local  maxima.  From  the  perspectives,  deviations  of  less  than  ten  degrees 
seem  to  be  veil  behaved  and  predictable  in  that  there  are  no  other  local  maxima  within 
a  10*  radius  of  any  given  maxima.  Convergence  limits  of  .001  degree  and  better  are 
realizable. 

3.  The  search  algorithm  does  not  provide  a  check  for  a  global  maximum.  This 
handicap  and  the  convergence  criterion  are  responsible  for  the  search  converging  at 
a  minima  adjacent  to  a  maxima.  This  tis  a  subject  of  interest  in  this  study. 

Results  on  a  test  case  are  presented  in  Appendix  E.  The  test  case  selected  is  for  motion 
in  an  orbit  inclined  at  43  degrees.  This  was  of  particular  interest  because  of  the  two 
adjacent  maxima  separated  by  a  minimum.  This  search  scheme  only  searches  for  the 
condition  in  which  the  first  derivative  is  equal  to  zero.  It  is  interesting  to  note  that  the 
figure  4 2D  (test  case.  i-43‘ ).  displays  17  different  areas  at  which  convergence  can  be 
achieved.  At  any  of  these  areas,  the  slopes  are  indeed  zero  and  the  search  scheme 


vould  these  out.  Achieving  convergence  does  not  signify  a  maximum  Aa.  It  is 

apparent  from  the  perspective  that  there  are  only  four  areas  in  vhich  Aa  is  maximum. 
Without  reference  to  the  perspective,  it  vould  be  most  difficult  to  relate  convergence 
values  to  maxima,  minima  or  inflection  points.  This  makes  it  so  imperative  that  the 
search  be  accomplished  hand-in-hand  vith  the  corresponding  perspective.  Only  then 
can  the  control  parameter  changes  be  meaningful. 


V.  Conclusion 
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Working  in  a  region  of  unknowns  and  searching  for  an  entity  yet  to  be 
described  is  just  what  was  initially  pursued  in  this  effort.  Given  a  multi-variable 

function  describing  the  change  in  the  semi-major  axis  (Aa)  with  an  established 
objective  of  finding  the  maximum  change  experienced  by  that  function  began  as  a 
formidable  task  of  high  interest  and  expectation.  The  pursuit  of  a  solution  using 
known  search  methods  e  g..  Stirling  Method,  to  evaluate  extrema  provided  much  hope 
for  convergence.  Futile  preliminary  attempts  to  identify  convergence  conditions  was 

due  to  poor  initial  guesses.  This  led  to  question  the  nature  of  the  Aa  function  at  hand. 
Understanding  what  the  function  does  is  quite  different  from  knowing  what  it  looks 

like.  The  big  question  is  this:  What  does  the  Aa  function  represent  graphically?  As  a 

resort  to  a  good  guess,  the  four-variable  Aa  function  was  graphed  perspective^  in 
three-dimensions  with  two  variables  held  fixed  and  two  varied  against  the  value  of  the 
function.  The  resulting  perspectives  represent  the  surface  of  the  function  for  the 
given  set  of "control”  conditions.  These  perspectives  give  a  good  overall  picture  of  the 

behavior  of  the  Aa  function.  Information  that  was  earlier  a  guess  can  now  be 
verified.  Regions  of  diminishing  returns  and  high  yield  can  be  identified  and  search 
patterns  can  be  concentrated  on  specific  areas.  Now  a  bound  exist.  Now  the  initial 
guess  can  be  nailed  down.  Search  schemes  which  are  sensitive  to  "good  initial  guesses" 
have  improved  reliability  for  convergence.  The  Stirling  Method  is  by  no  means 
excluded  here.  The  input  of  good  starting  points  (u‘)  to  initiate  the  search  for  the 

optimum  control  setting  that  yields  the  maximum  Aa  is  of  critical  importance. 

Convergence  in  two  or  three  iterations  shows  the  power  of  the  search  method  used. 
Moreover,  it  shows  the  goodness  of  the  starling  search  point.  Convergence  to 
coordinates  that  are  known  apriori  is  only  possible  from  reference  to  the  applicable 
perspectives  provided  that  the  starting  point  is  not  far  from  the  shown  perspective. 
Deviations  far  from  the  local  maximum  will  allow  the  search  mechanism  to  deviate 
away  from  the  intended  area  and  converge  on  another  unexpected  maximum.  This  is 
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normal  behavior  of  such  search  scheme  and  there  is  not  need  for  alarm.  The 
converged  control  settings  can  be  graphically  evaluated  by  the  use  of  the  perspectives 
defined  by  the  converged  control  values.  As  noted  earlier,  convergence  can  be 
achieved  at  any  area  satisfying  "zero  slope"  conditions  (e  g.,  marim*.  minima,  or 
inflection  points).  Only  by  looking  at  the  function's  perspectives  or  by  computing  its 
second  derivative  can  the  actual  maxima  desired  be  identified.  The  cross-relation 
between  the  search  scheme  results  and  the  perspectives  are  most  necessary  for  a 
meaningful  output. 

Comment  on  the  Equations  of  Motion 

Jenkins  [Ref:  101  provided  perturbation  equations  for  the  one-to-one  resonance 
case.  These  equations  have  been  thoroughly  verified  in  more  than  one  ways. 

1 .  The  equations  of  motions  were  integrated  for  the  resonant  case  (i.e.,  when 
the  orbital  mean  motion  is  equal  to  the  sail  precession  rate)  and  found  to  be 
correct. 

2.  When  the  equations  of  motion  were  used  to  generate  the  applicable  perspec¬ 
tives  at  a  90*  inclination,  it  was  discovered  that  at  a  coning  angle  0  -  35* .  a 
maximum  occurs.  This  is  in  direct  agreement  with  previous  results  arrived 
at  by  Tsu  [Ref:  21 1  and  Garwin  [Ref:  31  in  which  they  found  that  a  tilt  angle 
of  33  *  with  respect  to  the  sun  provides  optimum  change  in  the  semi-major 
axis.  The  constant  coning  angle  is  the  tilt  angle  of  the  sail  with  respect  to 
the  sun. 

These  indicate  that  the  perturbation  equations  arrived  at  by  Jenkins  do  indeed  describe 
the  appropriate  motion  of  a  freely  coning  solar  sail. 


Lessons  Learned 


It  is  fitting  that  whatever  was  learned  from  this  academic  effort  be  shared  with 
anyone  interested  in  pursuing  a  similar  quest.  The  work  done  here  is  without  its 
troubles  and  periods  of  despair.  The  biggest  difficulty  was  searching  for  the  maximum 
of  a  multi-variable  function  without  any  apriori  knowledge  of  its  behavior  under 
given  conditions.  The  suggestion  here  is  this:  represent  the  surface  (or  function)  in 
some  kind  of  perspective  (two-dimensional  or  three-dimensional)  and  observe  its 
behavior  as  certain  variables  are  changed  while  holding  others  fixed.  This  would 
provide  valuable  insight  into  what  can  be  expected  of  the  function.  Look  at  it  first  so 
that  they  are  no  surprises  later.  This  is  the  biggest  lesson  learned  from  this  academic 
effort. 


Recommend 


The  behavior  of  this  solar  sail  is  just  partially  known  by  the  exploitation  of  the 

Aa  function.  More  can  be  learned  by  numerically  looking  at  the  behavior  of  changes 
in  the  other  state  functions  eccentricity,  inclination  and  longitude  of  ascending  node. 
Although,  a  particular  aspect  of  the  inclination  is  addressed  here,  it  deserves  its  own 
segment.  The  numerical  evaluation  for  a  selected  performance  index  can  very  well  be 
supplemented  by  producing  the  three-dimensional  perspectives  associated  with  the 
case  of  interest. 

Another  pursuit  would  be  to  extend  the  single-stage  dynamic  system  to  a 
multi-stage  dynamic  system.  This  would  be  the  case  in  finding  the  optimum  series 
controls  which  would  generate  the  maximum  change  in  semi-major  axis  in  n 
revolutions.  This  would  entail  applying  the  single-stage  system  n-1  consecutive  times 
as  shown  previously  in  Figure  3  3  Such  a  multi-stage  system  can  have  practical 
applications. 
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Appendix  A 


Eclipsing  Effects 

The  eclipsing  or  shadowing  of  any  solar  radiation  dependent  spacecraft  is  an  important 
phenomenon  that  must  be  understood  and  compensated  for  in  the  design  stage  and/or 
by  some  control  mechanism  of  such  spacecraft.  For  this  solar  sail,  the  only  propulsion 
source  is  solar  radiation.  This  makes  shadowing  a  more  critical  phenomenon  for  a  solar 
sailcraft  than  for  a  spacecraft  with  a  variable- mass  propulsion  system.  To  achieve  the 
objective  in  maximising  select  orbital  parameters  (which  are  dependent  upon  the 
aaaount  of  solar  radiation  incident  on  the  surface),  the  aspects  of  shadowing.  e.g.. 

a.  when  does  shadowing  occur,  and 

b.  duration  of  shadowing. 

become  particularly  interesting  to  the  mission  designer  in  determining  the  optimum 
steering  controls  necessary.  The  dependence  of  the  changes  of  the  orbital  parameters 
on  the  shadow  time  for  this  particular  coning  solar  sail  can  be  readily  seen  in  the 
period  TP  term  in  Eqn  (221 ).  Note  that  TP  is  the  orbital  period  and  is  also  the  time  spent 
in  the  solar  radiation  environment  during  one  orbit.  Shadowing  would  result  in 
reducing  this  TP  value  and.  hence,  in  adjustments  to  the  amount  of  perturbations! 
changes  the  sail's  orbit  experiences.  Just  how  does  one  determine  a  and  b  above? 

Stoddard  (Ref:  19)  simplifies  the  aspects  of  shadowing  and  presents  a  means  of 
determining  whether  or  not  an  artificial  satellite  in  a  circular  orbit  is  shadowed  and.  if 
it  is.  what  its  duration  in  the  shadow.  The  "circular  cylindrical  shadow"  model  is 
employed  instead  of  a  "conical  shadow"  model  which  considers  the  umbra  and 
penumbra  shadow  components  (as  presented  by  Filler  in  Ref:  6)  to  arrive  at  the 
following  relationships  which  are  referenced  in  Figures  A1  and  A2: 


Fig.Al:  Earth's  Shadov  Geometry 


cos8=  ±  (l/e^)  11 -(Re/Rj)2]14 

(A- la) 

=  *  (seciju,,)  ll  -  (RE/R9)21 54 

(A- lb) 

0E-  2(180*  -  0) ; 

(A-2) 

(0b/36O*)TP; 

(A-3) 

where 

0£  *  geocentric  angle  of  travel  of  the  satellite  during  eclipse. 

isuN  9  geocentric  angle  between  the  sun  and  the  satellite's  orbit  plane. 

Note:  Front  Figure  Al.  this  is  equivalent  to  the  orbit  inclination  i. 

6  a  geocentric  angle  measured  in  the  orbit  plane  between  the  satellite  and  the 
conjunction  point  P  (Figure  A2). 

ejjj  a  eccentricity  of  the  elliptical  projection  of  the  earth's  shadow  on  the 
satellite's  orbit  plane. 

■•In:  This  can  be  geometrically  shown  to  be  equivalent  to  cos  i^yN.  This 
is  accomplished  in  (Ref:  19] 

Rjr  a  earth  radius  - 1  DU. 

R$  a  satellite  orbit  radius. 

TP  a  orbital  period. 

tgj  a  time  the  solar  sail  is  in  the  earth's  shadow  -  duration  of  shadow. 

The  limiting  case  between  eclipse  or  no  eclipse  is  obtained  from  Equation  (A-l )  when  0 
-180 degrees.  Therefore,... 


Rg  *  Rg esc isujj  -  Rgcscl 


(A-4) 


* 


Equation  (A-4)  is  valid  for  the  geometry  and  the  definition  of  the  inclination  of  the 
orbit  as  described  in  text.  That  is.  if  the  inclination  is  defined  as  the  inclination  of  the 
orbit  plane  wrt  the  ecliptic  plane,  then  this  Equation  (A-4)  "will  be  valid.  For  any 
application,  reference  to  the  Stoddard  article  is  highly  recommended. 

From  Equation  (A-4),  eclipsing  criteria  can  be  established  as  follows: 

if  Rj  <  R£csciSUN  -♦  a  shadow  will  take  place:  else  (A-3) 

if  Rj  >  REcsciSUN  =0  a  shadow  will  apt  take  place.  (A-6) 


A 

4 


4 


From  Equation  (A-3).  one  can  determine  that  any  artificial  satellite  in  the  ecliptic 
plane  (i  -  0* )  will  experience  shadowing.  Shadow  limits  for  various  inclination  i  and 
orbit  radius  R$  can  be  seen  from  the  plot  of  Equation  (A-3)  in  Figure  A3  above.  The 
area  to  the  right  of  the  "Limit  Line"  indicates  the  region  of  'NO  SHADOV. 


4 


A-4 


Rj  Cos  8 


0E*  TP(TU)  tjjjfTU)  ^(mim)  (%) 


changes.  From  Table  Al.  one  might  be  alarmed  at  the  increasing  magnitude  of  the 
shade*  time  tgj  for  increasing  orbit  radius  R$  and  its  resultant  decrease  in  these 

perturbations!  changes.  The  shadov  times  do  indeed  increase;  but.  the  percentage 
increase  of  shadov  time  over  the  period  TP  decreases  as  the  orbit  radius  increases. 
This  indicates  the  relative  impact  t^  has  for  high  orbiting  solar  sails.  The  further  out 

the  sailcrafl  is .  the  lesser  the  chances  of  shadoving  become.  For  the  solar  sailcraft 
initiating  its  maneuver  at  lov-earth-orbit  altitudes,  this  phenomenon  is  an  important 
issue. 

The  question  is  "ho*  much  shadoving  is  tolerable  vithout  considering  its  affects?' 
The  eclipse  factor  has  generally  been  used  to  ansver  this  question.  The  eclipse  factor 
is  defined  as  the  ratio  of  time  spent  in  the  sahdov  to  the  orbital  period  of  the  satellite. 
From  Eqn  (A-3).  this  can  be  eipressed  as 

Eclipse  Factor  -  t^/TP  »  0J/36O*  [Ref:  17]  (a-8) 

If  this  ratio  is  small,  then  shadoving  can  be  neglected.  Some  type  of  decision  criterion 
must  be  formulated  to  determine  an  ansver  to  this  question.  Such  criterion  will 
definitely  include  this  eclipse  factor  and  other  satellite/mission  dependent  parameters 
such  thermal  constraints,  etc. 

To  get  a  general  conceptual  feel  for  this  shadoving  effect  on  a  solar  sailcraft  in  the 
ecliptic  plane,  a  “simply  spinning”  case  vith  shadoving  and  another  vithout 
shadoving  vere  evaluated  and  contrasted.  The  data  in  Figure  A4  indicate  the  decrease 

in  total  change  in  ha  vhen  shadoving  effects  are  considered.  There  seems  to  be 
appreciable  relative  differences  betveen  these  computed  sets  of  values,  thus  indicating 
the  importance  of  the  shadov  phenomenon. 

Cases  dealing  in  the  ecliptic  plane  are  very  much  affected  by  this  phenomenon. 
Moving  out  in  altitude  and  to  inclinations  greater  than  zero  reduces  this  effect.  This 
study  avoids  this  phenomenon  by  having  the  initial  state  at  a  high  orbit  altitude  and  a 


non-ecliptic  inclination  such  that  condition  (A6)  is  satisfied. 

Escobal  and  Johnson  [Ref:  31  obtain  compact,  closed-form  expressions  for  the 
maximum  and  minimum  eclipse  durations  of  a  circular  orbit  with  knovn  semi-major 
axis  and  inclination  via  simple  geometric  constructions.  The  method  outline  provides 
the  designer  with  the  tool  to  evaluate  the  envelope  of  eclipse  durations  a  spacecraft 
would  experience  throughout  its  lifetime. 

Polyakhova  [Ref:  141  provides  a  more  extensive  and  thorough  coverage  of  the 
shadoving  phenomenon.  This  study  develops  the  solar  constant  equation  from  the 
basic  quantum  theory  of  light.  It  farther  investigates  the  shadow  effects  in  the  case  of 
the  radiation-pressure  influence  on  the  secular  acceleration  of  the  satellite,  i.e.,  on  the 

quantity  ATP/TP  (the  variation  of  the  satellite  period  during  one  orbit).  It  addresses 
orbits  of  arbitrary  eccentricity  and  develops  a  shadow  equation  and  provides  solution 
for  certain  simplifying  conditions.  The  author  develops  the  equation  using 
non-standard  reference  frames.  For  this  reason,  the  development  is  not  pursued  here. 
Strong  recommendations  is  made  to  it  for  those  pursuing  a  more  rigorous  approach  to 
this  shadow  phenomenon. 


fjp  mi  iwwiiivviwv';iirrv«’iu  .■  m  j  ■  j  ».■,■■.■  .vr?. 
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Appendii  B 


« 


Forming  tho  Sensitivity  Matrix  M 


O 


4 


% 


This  following  exercise  is  presented  to  bring  awareness  to  the  extend  of  labor  one  must 
suffer  in  determining  the  partials  of  the  change  in  the  semi-major  axis  da  to 

incremental  changes  in  the  control  parameters  (0,  Tf,  C,  PA).  Recall.  PA  s 
This  exercise  is  as  follows: 

IC]  *  [aAa/38  dda/dfi  dAa/ac  3  4a/dPA]  (B_,) 

where  (from  Equation  (2.21 ), 

da  =  a1-3  D  TP/n  [  dat  +  da2  4  da^  +  da^  +  da^  +  da6  ]  (B_2) 


Using  the  chain  rule  for  differentiation  and  looking  at  the  first  element  in  Equation 
(B-l). 


3  da 

a1-5  DTP' 

ddaj 
_  + 

dda2 

Ma3 
♦  _ 

3da4 
+  _  ♦ 

3da5 

adaj 
♦  _ 

39 

.  » 

.30 

as 

30 

30 

30 

30  . 

Working  solely  with  the  first  term  da^  on  the  rhs  of  Equation  (B-3), 

d«i  =  fc(dj)2  {(3^2  +  Pq) 009  ^  -  (3Pj  ♦  Pj)  sin  PA  }  .  (B-4) 


a 


B-1 


4 


adtj/ae  -  d(a{)2m  ((3p2  +  cos  pa  -  (3^  ♦  p5)  sin  pa  } 

+  (d|  )2  {3ap2/a0  ♦  apq/ao }  cos  pa 
-  (d^2  [3d$xm  +ap5/a0}  sinPA 

Expanding  the  d  j  and  d§'s  terms  above,  where... 

d]  -  sin  8  cos  <  cos  i, 

and 

Pi  »  sin  8  cos  C. 

P2  *  sin  0  sin  <  cos  i\, 

P4  =  sin  0  sin  C. 

-  sin  8  cos  ti, 

one  arrives  at  the  following  results: 


(B-5) 


(B-6) 


(B-7a) 

(B-7b) 

(B-7c) 

(B-7d) 


ada,/a0  -  1/8  [  2s20  c2C  A  ((3s0s£ct)  +  s8s<)cPA  -  (3s0c<  -  s8cti)sPa} 
♦  ( 1  -  c20)  {(3c0s<cn  ♦  c8sC)cPA  -  (3c0c<  ♦  c0cti)sPa}  ]  (B-8) 


where  c  a  cos  and  S  a  sin. 


Note  that  Equation  (B-8)  takes  care  of  the  first  term  in  Equation  (B-3).  There  are 
several  more  terms  that  require  expansions  and  differentiations  wrt  the  control 

parameter  0.  Specifically. 


dd&,/d0  =  ? 


dda*/dO  -  ? 


Extending  this  to  the  partiais  of  the  remaining  control  parameters  (71,  (.  and  PA),  one 
can  easily  imagine  the  amount  of  labor  required  and  the  countless  room  for  error(s). 


dA3  /dt]  *  ?  dAa  /d£  *  2  3A a  /dPA  *  2 


This  study  circumvents  this  potentially  troublesome  area  and  approaches  the 
determination  of  the  Sensitivity  Matrix  C  by  employing  numerical  techniques.  This 
approach  is  discussed  in  the  text. 


The  Surface  Program.  The  Fortran  Program  constructed  to  generate  the  various 
3-dimensional  perspectives  of  4a  is  included  here  in  its  entirety.  It  is  designed  to  be 
user  friendly;  it  Till  lead  any  prospective  user  to  easily  generate  similar  perspectives 
provided  that  there  is  some  basic  understanding  of  the  mechanics  of  the  UVIX/VAX 
systems  and  the  operation  of  the  HP7220v  Plotting  Table. 

Generating  the  Data.  Equations  (221)  ,  (222),  (229)  and  (2.30)  vere 
programmed  in  a  user-interactive  program  called  ZDATA  to  generate  the  changes  in 
the  semi-major  axis  parameter.  Desired  inputs  queried  include  the  tailoring: 

a.  the  initial  semi-major  axis 

b.  the  orbit  inclination 

c.  the  angular  momentum  orientation  angles  ( r?  and  O 

d.  the  range  in  the  coning  angle  (0),  and 

e.  the  range  in  the  phase  difference  angle  (PA  *  <£-(/'). 

Additional  inputs  in  scale  factor  are  required  to  arrive  at  a  data  value  that  can  be  easily 
interpreted.  Scaler  increments  the  ranges  of  the  angles  0  and  PA  to  allovr  a 
macro-perspective  or  a  micro-perspective  of  the  a  function.  Scales  raises  the  data  to 
a  value  between  -IQjQ  and  10j0  tar  easy  interpretation  of  relative  magnitudes  of  4a. 
The  actual  4a  is  preserved  in  the  parameter  called  z. 

The  information  generated  is  placed  into  a  data  file  (called  plotdata)  compatible  to  the 
plotting  routine.  The  Program  Listing  is  attached. 

The  Plotting  Routine.  This  is  basically  the  set  of  instructions  tar  getting  the 
data  plotted  on  an  HP7220V  Plotting  Table.  Additional  information  can  be  obtained  from 
the  S-Package  routine.  The  instructions  are  as  tailors: 


1.  Locate  the  terminal  co-located  Tith  the  HP7220v  Ploting  Table. 

2.  Login  as  usual  and  move  into  the  directory  in  Thich  the  plot  data  is  located. 

3.  On  the  tenninaL  enter  the  toUoring: 

a.  set  term=hl9  (this  connects  terminal  to  UH/VAX  system) 

b.  bit  SETUP  key. 

c.  hit  B  key  and  set  baud  rate  to  2100;  'verify  similar  setting  on  rear 
of  HP7220v  device. 

d.  hit  SETUP  key  to  lock  value  in. 

1.  Turn  plotter  on  and  set  paper  &  pens  up;  place  plotter  on  STBT 

3.  On  the  terminal  enter  the  toUoring  (note:  %  and  <  are  computer  prompts) 

a.  %S  (response Till  be <) 

b.  <  hp7220v  (this  identifies  plotting  device) 

c.  >  i-matrix(read(“daia  fi1ename*)^mre*xx.bynrw*T) 

There  eiMmnm  »  WtoBMM  imM  to  pnnwwHtig  am 

nror*  ranges  scales  +1 

(see  S-Package  rafernca  tor  additional  options) 

d.  z_s-min(s) 

e.  z_z/max(z) 

f.  perspCt)  (see  S-Package  reference  tor  additional  options) 

Mote:  3d  and  3e  dimension  the  data  to  the  maximum  value  of  x  contained  in  data  set. 

Relative  magnitudes  can  be  expressed  by  inserting  the  maximum  value  of  4a 
obtainable  in  step  3d  above.  The  resulting  perspective  Till  give  the  relative 

magnitudes  of  the  a  values  art  to  this  maximum. 

Caution:  The  underscore  _  tolloring  the  term  x  in  the  instruction  above  is  an 
important  function  of  the  plotting  routine  and  must  not  be  neglected. 


Proqron  Listing 


consent:  This  progran  generates  the  plotdata  for  input  into  the 
c  S-Package  routine.  The  data  is  used  in  progran  persp(z). 
c  eodified  for  pa  us  theta  (pa  is  horizontal) 
c  surface  generation  routine  (2  uar table;  3-die) 
c  Double  precision  is  not  required  for  perspectiues. 
c 

integer  i, j ,np1,np2,scalex 

integer  dif1,dif2, i 1 ,  i2,f1,f2,dect, iter 

i nt eger  dec3 , decl , i case 

real  da(0: 100,0: 100)jdax(0: 100,0: 100),scalez,datot,coeff 
real  a, ix, in,tp,sac,conu,pi 
real  theta, eta, etax, chij chi x, pa 
real  b1,b2,b3,b1,b5,b6,b?,b8,d1,d2,d3 
real  da1,da2,da3,da1,da5,da6 
c 

data  sac/. 00000165/ 


consent :  open  up  tapes  file  for  output.  Plot  data  contains  the 
c  all  the  necesarry  data  to  generate  the  perspectives, 

c  Ho  consents  are  alloned  in  this  data  set.  Plot  info  con- 

c  tains  all  pertinent  infornation  to  plotdata. 


open(un i t-11 , f i I e- ' p I ot dot a ’ , access" ' sequent i a I ' , 
$  st at us"' nee') 

open(unit"12, f i leB< plot  info', access" 'sequential ', 
$  status"'nen>) 
renind(unit*l 1) 
renind(unit"!2) 


10  print*, 'Enter  a, i, eta, chi  (canonical  &  degrees)  (real)' 
read*, a, ix,etax,chix 

print*,  'CAUTION:  Do  not  exceed  a  100  X  100  array!' 
print*, ‘Enter  SGBLEX  (e.g.,  5  pt*deg  ■  5)  (integer)' 
read*, sea I  ex 

20  print*, 'Enter  IHITIOL  &  FINRL  points  for  THETR  (integer)' 
read*, i 1 , f 1 

print*, 'Enter  IHITIOL  fc  FIHflL  points  for  PO  (integer)' 
read*, i2,f2 

print*, ‘Enter  SCflLEZ  (desired  height  for  OB)  (real)' 
print*, 'note:  scalez  »  10000.  (das  right!  10  K)  eorks' 
read*, sea I ez 

print*,  'Bre  you  sure,  Solar  Sailor?  1**YES;  0-N0' 

read*, dec 1 

if(dec1 .ne.1)  go  to  10 

c _ 

coeeent:  Coepute  soee  ieportant  paraeeters. 
c 

pi  -  1.*atan(1 . ) 
tp  -  2.*pi*sqrt(a*a*a) 
cony  ■  pi/180, 
in  ■  ix*conv 
chi  ■  chix*conv 
eta  ■  etax*conv 
c 

coeeent:  echo  input  back  to  screen. 

print*, 'period  tp  (tu)  "^tp 

print*, ' incl inat ion  (deg)  a',ix 

print*, 'eta  (deg)  “^etox 

print*, 'chi  (deg)  "',chlx 


c 


k 

L 


o 


I 

I 


4 


4 


coment:  Determine  the  intervals  for  plot.  Intervals  eust  be 
c  equally  spaced, 

c 

difl  -  Iabs(f1  -II) 
dif2  -  iabs(f2  -  12) 
npl  -  difl/scalex 
np2  »  dif2/scalex 
c 

print*, 'difl  -  ',dlf1, *  dif2  -  ',dif2 

print*, 'npl  -  ',np1,'  np2  *  *,np2 

iter  -  0 
c 

•rite(12,190)a, ix,etax,chix, i 1 , f 1 , i2,f2,scalex,scalez 
c 

•rite(!2,l80) 

•rite(*,180) 

180  fornat (//30x, ' Pfl ' , /) 

c _ 

c**************************  oo  LOOP  ************************** 

c 

c  npl  is  the  •  of  points  in  Theta  paraneter. 

c  np2  is  the  •  of  points  in  Pfl  paraneter. 

c  scalex  is  the  9  of  divisions  desired, 

c  scalez  is  the  factor  required  to  raise  the  0a  value  to  sone 

c  value  that  can  be  eorked  eith. 

c 

do  100  i  ■  0,np1 
theta  *  (il  ♦  i*scalex)*conv 
c  print*,'  theta  -  ', theta 
do  110  j  •  0,np2 
pa-  (i2  -  j*scalex)*conv 
c  print*,'  pa  ■  ',pa 

c 


4 


C-5 


coaaent :  Coapute  the  various  parameters  in  the  Da  equation. 

bl  •  sin(theta)*cos(chi) 
b2  •  sin(theta)*sin(chi)*cos(eta) 
b3  *  sin(theta)*sin(chi)*cos(eta) 
bl  *  sin(theta)*sin(chi ) 
b5  *  sin(theta)*cos(eta) 
b6  •  cos(theta)*cos(chi )*sin(eta) 
b?  ■  sin(theta)*sin(eta) 
b8  *  cos(theta)*cos(eta) 

dl  -  sin(theta)*cos(chi )*cos( in) 

d2  ■  sin(theta)*sin(chi )*cos(eta)*cos( in) 

♦  sin(theta)*sin(eta)*sin( in) 

d3  -  cos(theta)*sin(chi)*sin(eta)*cos(in) 

-  cos(theta)*cos(eta)*sin(in) 

dal  -  . 25*d1 **2 . *( (b2*3 . *b1)*cos(pa)-(3 . *b1 ♦b5)*s i n(pa) ) 

da2  -  - . 50*d1 *d2*( (b2-b1)*si n(pa)+(b5-bi )*cos(pa) ) 

da3  -  2 . 00*dl *d3*(b6*cos(pa)-b3*s i n(pa) ) 

dal  •  2.00*d2*d3*(b3*cos(pa)+b6*sin(pa)) 

daS  -  . 25*d2**2 . *( (3 . *b2+M )*cos(pa)-(b1 *3 . *b5)*s i n(pa) ) 

da6  ■  d3**2.*((b2+b4)*cos(pa)-(b1+b5)*sin(po)) 

datot  *  da1*da2+da3*dai+da5*da6 
coeff  •  a**1 .5*sac*tp 
dax(i,j)  -  coeff*datot 
da(i,j)  -  scalez*dax(i ,j) 
iter  *  iter  ♦  1 
110  continue 


erite(1 1 ,200)(da( i , j ), j“0,np2) 
•rite(12f200)(dQ(iJ),j-0,np2) 

•rite(*,200)(da(i, j), j-0,np2) 
tOO  continue 

_ FORftATS _ 

170  foreat</5x,'tp  •' ,ib.Zt2x,  '  in  ’eta  -,Jf6.3/chi 

$  f6.3) 

% 

190  for«at(//, 'a  ■' ,  f6.3,2x, '  in  *' ,  f6.3,2x,'eta  ■',f6.3,2x,  *chi 
$*',f6.3,  2x, ‘theta  froes'dl,'  to' ,  i4,2x,  'pa  froe: 

$  '  1 14,  *  to',  i4, /,  'seal  ex  -',f3.l!2x,  'scalez  -  * ,  f  1 0 . 2) 


200  foriat(100(f7.3)) 


print*, ’Do  you  eant  to  try  another  caee????  1mYES;  0-N0 
read*! icase 
i f ( icase.ne.O)  then 
print*,  'a  “',a 
print*,' i  -',ix 
print*,  'eta  -',etax 
print*, 'chi  "'jChix 

print*, 'Any  change  in  a, i, eta, chi?  1»VES;  0-N1 
read*,dec3 

i f(dec3.ne.O)  go  to  10 

print*, 'theta  starts  at ' , i 1 , *  and  ends  at',f1 
print*, 'pa  starts  at  * , i 2, '  and  ends  at',f2 
print*, 'Any  change  in  theta  and  pa?  1*YES;  0- 
read*, decl 

i f(decl.ne.O)  go  to  20 

endi  f 


endfi le(unit-ll) 
endf i le(unit*12) 
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Appendix  D 
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Contents.  The  tolhrire  constructed  to  determine  the  set  of  control  vectors  thst 
would  make  the  performance  index  J  stationary  consisted  of  the  foiloxring. 


1.  Main  Program:  ZERO  f 

2.  Subroutine: 

ECHO 

3.  Subroutine: 

UPART 

4.  Subroutine: 

DELTA 

3.  Subroutine: 

INV4X4 

Description.  Here  is  a  short  synopsis  of  the  function  of  each  program  used  in 
this  search  scheme.  The  entire  listing  is  documented  further  in  each  section  of  the 
■•in  where  additional  information  would  facilitate  its  understanding.  The  programs 
are  as  follows: 

Pmtftin  7EROF  This  is  the  main  driving  program  that  employs  the 
above  subroutines  to  develop  the  necessary  first  and  second  partials  and  to  determine 
the  zero  of  the  optimality  condition.  It  is  a  user  friendly  program  written  to  work 
interactively.  It  queries  the  user  for  necessary  input  data,  an  iteration  limit,  and  a 
convergence  limit. 

Subroutine  ECHO.  The  input  data  along  with  generated  preliminary 
information  (e.g.  orbital  period)  are  echoed  back  to  the  user  for  verification. 


Subroutine  UP  ART  The  partials  of  the  state  function  Fa  vith  respect  to 

the  control  variables  (8.1).  (.  $-\R)  are  computed  via  the  Stirling  (a.Jt.a  modified 
Nevton-Rhapson  or  central  differencing)  iteration  method.  This  program  is 
summoned  four  times  to  construct  the  optimality  control  vector  and  supporting 

partials  (PC]  Output  is  d%t/(dUjUj^  in  the  form  (PU A]  and  (PUB). 

Subroutine  DELTA.  This  uses  the  perturbation  equations  to  determine  the 
changes  in  the  orbital  parameters  a.  Output  is  Aa. 

Subroutine  INV4I4.  This  is  used  to  invert  the  4X4  (PC)  matrix  so  that  the 
8u  required  can  be  determined.  It  outputs  the  determinant  and  inverse  of  (PC)  -  (PCI). 

Program  Listing  Attached  is  the  entire  program  listing  of  the  main  program 
ZERO  J.  All  supporting  subroutines  required  to  compile  and  execute  the  main  program 
are  included  here;  it  is  a  stand-alone  program. 


The  foil 00 i ng  are  the  computer  codes  eritten  to  search  for  the 
stationary  value  of  the  Performance  Index  J.  The  Stirling  He t hod 
is  used  to  iterate  and  find  the  zero  of  the  function  that  eould 
alloe  this  stationary  value  to  exist.  Equations  found  in  Section  3 
under  DISCUSSION  are  the  basis  for  this  program. 


PR06RRR  ZERO 

C 

C  Last  modified  on  25  Nov  84  0600  hrs 

C 

C  This  program  drives  the  partial  routines  to  determine  the  zero 
C  of  the  optimality  condition: 

C  Lamda  ( Transpose ) *Upart ia Is  -  [0] 

C 

INTEGER  DECO, 0EC1 , DEC3, 0EC5, DEC6, DEC?, DECLRS 
I NTEGER  DEC8, 0EC1 0, DECI 3, DEC1 4, DEC50, DATE 
INTEGER  iterx, i i ,  j j ,m,n,part,k, iter, case 
C 

DOUBLE  PRECISION  PU(4,4),PC(4,4),PCI<4,4),DIFF(4,4) 

DOUBLE  PRECISION  del (4),deldeg(4),C(4),PUA(4,4),PUB(4,4) 
DOUBLE  PRECISION  a, i ,ax, ix,cv,det,tp,sac,mu,dap,dam,da 
DOUBLE  PRECISION  theta, eta, chi, pa, thetax,etax,chix, pax 
DOUBLE  PRECISION  conv,dconv,tconv,pi ,tpi 
DOUBLE  PRECISION  dthe,deta,dchi ,dpa,delu,du,dux 

DATA  ax , i x/2 . d+00 , 30 . d+00/ 

DATA  t  het  ax , et  ax , ch i x , pax/90 . d+00 , 0 . d+00 , 0 . d+00 , -90 . d+00/ 
DATA  SRC, i ter/4. 65d-06,0/ 

DATA  0C0HU/7 . 9053682d+00/ 

ORTfl  TC0HU/806 .Bill 8742d+00/ 

DATA  mu, cv/1 .d+00, .001 d+00/ 

0ATA  dap, dam/0. d+00, 0. d+00/ 

DATA  dt he, det a/0. d+00, 0. d+00/ 

ORTA  dchi,dpa/0. d+00, 0. d+00/ 

DATA  de I u, dux/. 01 d+00, .05d+00/ 


COntlEHT:  Open  up  a  tapes  for  print  files  to  be  stored  in. 

C  unit  19  contains  the  del(ii)  results  eith  sys  description. 

open(un i t-1 9, f i I e- ' output  19 1 , access- ' sequent i a  I ' , 

$  status-' nee 1 ) 
reeind(unit-19) 


COHHENT:  This  starts  the  interact iue  eode  for  data  input  purposes 
C 

print*/ Enter  DATE  and  CASE  HO  (e.g.  180984,5):' 
read*, date, case 

204  print*, 'Any  changes  in  OABITAL  PAAAHETEAS?  1-V;  0-H' 

read*,  DECO 
if(DECO.ne.O)  then 

print*, 'See i-eaj or  Axis  a  -',ax  /Any  Change?  1-V;  OH 
read*, DEC 1 
i f(0EC1 .ne.O)  then 

print*,'ENTEB  a  (Oil)  (double  precision)' 
read*, ax 
end  if 

print*/  Inc  I  inat  ion  i  -',ix  /Any  Change?  1-V;  OH 

read*,DEC3 
i f(OEC3.ne.O)  then 

print*,'ENTEA  i  (deg)  (double  precision)' 
read*, ix 
endi  f 
endi  f 


208  cont i nue 
C _ 

print*, 'Current  theta  -',thetax/  AHV  CHANGE?  1-V;  0-H 
read*,DEC5 
i f(DEC5.ne.O)  then 

print*/Enter  theta  (deg)  (double  precision)' 
read*,thetax 
endi  f 

print* /Current  eta  -',etax,' 
read*,  0EC6 


AHV  CHANGE?  1-V; 0-H' 


i f(DEC6.ne.O)  then 

print*/ Enter  eta  (deg)  (double  precision)' 

read*, et ax 
end  if 

print* /Current  chi  -',chix/  flNV  CHANGE?  1-V;0-H' 

read*,  DEC7 
i f(DEC7.ne.O)  then 

print*/ Enter  chi  (deg)  (double  precision)' 

read*, chi x 
end  i  f 

pr int*/ Current  pa  -', pax,'  flHV  CHANGE?  1-VjO-M* 

read*,  DEC8 
i f(OEC8.ne.O)  then 

print*/ Enter  pa  (deg)  (double  precision)' 
read*, pax 
end  if 

print*,' Current  de I u  - ' , de I u/  RHV  CHANGE?  1 -V  j 0-H ' 

read*, dec 10 
if(declO.ne.O)  then 

Print*/Enter  delu  (double  precision)' 

read*,delu 
end  if 

print*/ Current  du  -',dux/  RHV  CHANGE?  1-V;  0-H' 

read*,dec13 
i f(dec13.ne.O)  then 

print*, 'Enter  du  (double  precision)' 
read*, dux 
endi  f 

print*,'Do  you  eant  to  iterate?  1-V;  0-H* 

read*,decM 
if(decH.ne.O)  then 

print*,'Enter  Ho.  of  Iteration  (integer  fonat)’ 
read*, iterx 

print* /Convergence  Lieit  -',cv,'  RHV  CHRHGE?  1«V;,0-N’ 
read*,dec50 

i f(dec50.ne.0)  then 

print*,'Enter  Convergence  Lieit  v  (double  precision)' 
read*,cv 


COffflENT:  Initialize  the  systea  state  paraaeter: 

a  ■  ax 
C 

C0M1EHT :  CRLCULRTE  PI  IH  R  UNIQUE  URV. 

C 

PI  -  4.d+00*DRTRN(1.d+00) 

TP  I  -  2 . d+00*P I 
COHU  -  PI/180.d+00 
i  *  ix*conv 
C 

C0I1I1EHT:  Calc  the  period  of  the  initial  systea  a/o  perturbation. 

C 

TP  -  TP  I *DSQRT (R**3 . d+00/HU**2 . ) 

C _ 

209  cont i nue 
c _ 

coaaent :  If  iteration  is  desired,  209  ail  I  bring  in  the  Unea  values 
c  that  aere  coaputed  froa  UPRRT.  The  control  paraaeters  are 

c  updated  so  that  they  are  stepped  through  the  function 

c  until  convergence  is  aet  or  until  an  iteration  I  in  it 

c  i s  reached .  Reca 1 1 , 

c  Oelta  u  ■  Unea  -  Uold 

c  Unea  *  Uold  +  Delat  u 


COnnEHT:  PRINT  HERDER  FOR  OUTPUT; 


i f( I .EQ.0.d+00)  then 
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CALL  ECHO (fix, lx,TP,etax,chix,thetax,pax) 


COflUEHT :  COHUERT  RLL  DEGREES  TO  RRDIRNS: 

C 

theta  -  thetax*conu 
eta  ■  etax*conu 
chi  ■  chix*conv 
pa  *  pax*conu 
C 

COnflEHT:  This  1st  call  to  UPRRT  calculates  the  C  Hatrix. 

C  ****H0  INCREHENTAL  CHARGES  ARE  INTRODUCED  HERE**** 


CALL  UPART(a, i , theta, eta, chi ,pa,tp,tpi ,conu,c, 

$  de  I  u,  dt  he,  det  a,  dch  i ,  dpa,  du,  pu,  eu,  det ,  sac,  ■,  n,  part ) 


COffflENT:  The  required  aatrix  ie  [partial  F/partial  u]. 

C  This  aatrix  is  constructed  a  roe  at  a  tiae. 

C  Any  subsequent  call  to  UPRRT  is  to  build  the  ffatrix  PC 

C  eleaents  by  roes.  This  is  very  IHPORTRHT.  UPRRT  is 

C  called  taice  aore  to  fill  this  required  PC  aatrix. 

C 

C  part  -  2  denotes  the  2nd  partials  are  being  coeputed. 

c 

coaaent:  Set  up  the  initial  PU  aatrix  and  populate  it  eith  zeros, 
c  This  aatrix  is  used  later  but  filled  in  by  roes  at  a  tiae. 

c  This  population  of  zeros  is  required  to  prevent 

c  segmentation  errors.  recoil:  [PU]  »  [1x4]. 

c 

do  90  i i  ■  1,4 

do  92  jj  -  1,1 
pud  i,jj)  ■  0.d+00 
92  continue 

90  cont i nue 
c 

coaaent:  This  next  section  contains  the  necessary  do  loops  to  build 
the  PUR  abnd  PUB  aatrices.  THe  loop  continues  in  the 
next  feo  pages. 


c 

c 


© 


c 


© 


4 


* 


Q*************************  qo  LOOP  ****************4*************** 
conent :  the  k  loop  is  responsible  for  computing  PUR  and  PUB. 
c 

part  *  2 
c 

do  1000  k  ■  1,2,1 

I f (k . gt . 2)  go  to  1000 
if(k.eq.l)  du  -  dux 
if(k.eq.2)  du  a  -dux 
c 

conent:  the  ■  loop  is  responsible  for  generating  the  eleeents. 
c 

do  2000  e  ■  1,4 
c 

if(e.eq.l)  then 

theta  *  thetax*conv  ♦  du 
eta  a  etax*conv 

chi  ■  chix*conv 

pa  a  pax*conv 

else! f(e.eq.2)  then 
theta  ■  thetax*conv 
eta  ■  etax*conu  ♦  du 
chi  a  chix*conv 

pa  *  pax*conu 

elsei f(e.eq.3)  then 
theta  ■  thetax*conu 
eta  a  etax*conv 
chi  ■  chix*conu  ♦  du 
pa  *  pax*conv 
elsei f(e.eq.l)  then 
theta  ■  thetax*conv 
eta  a  etax*conw 
chi  ■  ch i x*conv 
pa  ■  pax*conv  ♦  du 
end  i  f 


i 
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COtffiENT:  This  next  call  to  UPIIRT  in  the  LOOP  generates  the 
elements. 

C  if  k  ■  1  then  PUR  is  being  filled. 

C  if  k  ■  2  then  PUB  is  being  filled. 

r 


cal  I  UPRRT(a,i  , theta, eta, chi , pa, tp,tpi ,conu,c, 

$  de I u, dt he, deta, dch i , dpa, du , pu,  mi, det , sac, e, nf part ) 


cont i nue 


if(k.eq.l)  then 

or ite(1 8,562)  k 
orite(19,562)  k 
do  40  ii  -  1,4 

do  42  jj  -  1,4 
pua(ii,jj)  -  pu(ii,jj) 
continue 

cont i nue 

erite(18,910)((pua(i i,jj),jj-1,4), i i-1,4) 
orite(19,910)((pua(i i,jj),jj-1,4), i i-1,4) 

elsei f(k.eq.2)  then 

erite(!8,564)  k 
orite(19,564)  k 
do  44  ii  -  1,4 

do  46  jj  -  1,4 

pub< i i , j j )  -  pu ( i i , j  j ) 

cont i nue 

cont inue 

orite(18,910)((pub(ii,jj),jj-1,4),i i-1,4) 
erite(19,910)((pub(ii,jj),jj«1,4),i i-1,4) 

end  if 
continue 
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COfWEHT;  Hoe,  get  the  difference  [PUR  -  PUB};  divide  by  2  x  du. 
C 

erite(19,566)  k 
c 

do  21  ii  -  1,4 
do  22  jj  -  1,1 

diff(ii,jj)  ■  puo(ii,jj)  -  pub(ii,jj) 

22  cont i nue 

21  continue 
c 

erite(19,910)((diff<ii,jj),jj-1,1),i i-1,1) 
c 

c  _  Rsseeble  the  PC  eatrix. _ 

c 

erite(19,568)  k 
c 

do  25  ii  -  1,1 
do  26  jj  -  1,1 

pc( i i , j j )  -  diff(ii,jj)/(2.d*00*du) 

26  cont i nue 

25  cont i nue 
c 

erite(19,910)((pc(i i, j j), j j-1,1), i i-1,1) 
c 

c  Calculate  the  INUERSE  of  eatrix  PC 

c _ 

call  IH1Mx4(pc,pci ,det) 


■r i te( 19,912)  det 
c 

coeeent:  Tereinate  Job  if  eatrix  PC  is  singular, 
c 

i f(det .eq.0.0d+00)  go  to  999 
c 

erite(19,580) 

erite(19,910)((pci(ii,jj),jj-1,1),i i-1,1) 
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COIfflENT:  Determine  the  delta  u  requried  and  check  for  convergence. 
C  If  no  convergence,  repeat  search  for  a  zero  by  altering  a 

C  selected  control  paraaeter. 

C 

■rite(*,590) 

arite(19,590) 

C 

COIfflENT:  Intiolize  the  DEL  aatrix  to  zero:  [1x4]  *  [0] 

C 

do  800  ii  -  1,4 

del ( i i)  -  0.0d+00 
800  cont i nue 


C 

C0HHENT:  Coapute  delta(u):  delta  u  •  U(nea)  -  U(old) 

C  delta  u  -  del(i)  -  -  C(1x4)*PCI[4x1] 

C  Unee  -  Uold  ♦  delta  u 

C 

COIfflENT:  If  deltau  is  not  approx  0.0,  then  try  it  again  eith  a 
C  a  slightly  different  choice  of  orientation  angle. 

C  Add  the  appropriate  delta  u  to  strat  the  iteration. 

C 

do  801  i i  *  1 ,4 
do  802  jj  -  1,4 

del(ii)  -  -(C(jj)*pci(jj, i i)  ♦  del(ii)) 
deldeg(ii)  ■  del(ii)/conv 
c 

802  cont i nue 

801  continue 


c 


coaaent:  Print  out  the  required  change  in  u  in  degrees  and  radians. 


6 


•rite<*,590) 

erite(19,590) 

erite(*,910)  del(1),del(2),del(3),del(1) 
erite(19,910)  del (I ),del(2),del(3),del<1) 
erite(*,592) 

•r I te(1 9, 592) 

•rite(*,910)  deldegd  ),deldeg(2),deldeg(3),deldeg(4) 
erite(19,910)  deldegd ), del deg(2), del deg(3), del deg(1) 

c _ 

coeeent :  Detereine  if  convergence  is  set. 
c 

c  if(del(l). It. cv.and.de  I (2). It .cv. and. del (3). It.cv. 

$  and. del (1). It.cv)  then 
c 

•rite(*,990) 
erite(19J990) 
end  if 

C _ 

COnnEHT;  Set  up  the  iteration  block: 

C 

i fdter.eq.  iterx)  then 


!• 

theta 

•  theta  ♦  del(1) 

eta 

■  eta  ♦  del (2) 

chi 

-  chi  ♦  del (3) 

pa 

■pa  ♦  del(4) 

% 

C 

thetax 

■  thetax  ♦  deldeg(l) 

etax 

■  etax  ♦  deldeg(2) 

chix 

■  chix  ♦  deldeg(3) 

pax 

■  pax  ♦  deldeg(l) 

• 

C 

cal  1  echo (ax, ix,tp, etax, chix, thetax, pax) 

cal  1  delta(a, i ,eu, theta, eta, chi ,pa, tpi ,sac,tp,da) 

c 

erite(l9,915)  da 

erite( 

*,915)  da 
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go  to  999 

elsei f(decll.ne.O)  then 
iter  ■  iter  ♦  1 
c 

coeeent :  echo  the  results  of  the  run. 


cal  I  echo (ax, ix,tp,etax,chix,thetax,pax) 


go  to  209 


endi  f 

999  cont i nue 


endf i I e(un i t—1 9) 


FORURTS 


100  FORflRT  ( 1  HI ,  /5X, '  HOT  I  OH  IH  THE  ECLIPTIC  PLRNEMX, 

S’DRTE: ', i8,3x, ’CASE; *, i3,3x, ’ ITEB; iO,/) 

101  F0RnflT(1H1,/5X/M0TI0N  IN  THE  HON-ECLIPTIC  PLRHE\5X, 

$ ’ DATE : ' , i6,3x, 'CASE: * , i3,3x, *  ITER: * i3,/) 

515  foreat (/39x, ‘Uector  C  [1x1]') 

562  foreat (/30x, 1  For  k  ■  ' ,  U  /  #  Matrix  PUR  [1x1]',/) 

561  foreat (/30x, ‘For  k  ■  ',11,',  Matrix  PUB  [1x1]*,/) 

566  foreat  (/39X|  'Matrix  DIFF  [1x1]  ‘,20x/k  -  Ml,/) 

568  foreat (/39x, ' Matrix  PC  [1x1]  ‘ ,20x, ‘k  -  ‘ , il,/) 

580  foreat (/39x, 1  Matrix  PCI  [1x1]', 20x, ‘k  -  * , il,/) 

590  foreat (/, 'Uector  DEL  (rads)  [1x1]  ’,/) 

592  foreat  (//Uector  DEL  (degs)  [1x1]  ’,/) 

910  foreat (1(d20. 10)) 

912  foreat (/5x,'0et ere inant  ■  ',f50.10) 

911  foreat (/5x, f20. 10, f20. 10, f20. 10, f20. 10) 

915  foreat (/5x, '0a  ■  ',f20.10) 

990  foreat (//5x, 'By  joue,  I  think  you  got  it,  Solar  Sailor1) 


SUBROUTINE  UPRRKa,  ijthetojetajchijpajtp^pijconv^c., 

$  delu,dthe,deta, dchi , dpa,du,pu,au, det, sac, «,n, part ) 


c 

c  This  subroutine  calculates  tlatrix  PU.  It  uses  the  DELTA 

c  routine  to  calculate  the  changes  in  the  orbital 

c  parameters, 

c 

integer  a,  n,  part 

double  precision  a, i ,sac,tp,tpi ,au,conv,di f 
double  precision  thetajetajChi ,pa,dthe,deta,dchi ,dpa 
double  precision  thetapjetapjChipjpapjdelujdu 
double  precision  thetaa^etaa, chin, pan 
double  precision  detjdapjdaa 
double  precision  PU(1, 4),C(4) 

c - 

coaaent:  Initialize  the  da  value  to  zero, 
c 

dap  a  0.d+00 
doa  *  0.d+00 
c 

coaaent:  Get  the  orbital  parameter  perturbat ions  for  the  given  set 
c  of  orbital  parameters  and  orientation  angles.  Use  these  for 

c  the  basis, 

c 

coaaent:  This  first  call  to  DELTR  coaputes  the  F(old)  values, 
c  Ue  aant  to  deteraine  the  perturbations  eithout  small 

c  changes  dt he, det a , etc.  i.e.,  dthe-deta-dchi-dpa-O.O 

c 

coaaent:  The  next  page  contains  the  00  LOOP  for  the  construction  of 
c  components  of  the  PUfl  and  PUB  matrices.  This  constitutes 

c  the  n  loop  of  the  (a,n)  LOOP  started  in  the  HR  I H  program, 

c 
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DO  LOOP 

coeeent;  Determine  the  part  idle.  Start  a  do  loop  for  the  eth  and 
c  nth  components  of  the  eatrix.  This  is  done  by  roes, 

c 

i f(part .eq. 1 )  ■  -  1 
do  400  n  -  1,4 
c 

if(n.eq.l)  then 
dthe  ■  delu 
deta  -  0.d+00 
dch i  -  0 . d+00 
dpa  -  0 . d+00 
elsei f(n.eq.2)  then 
dthe  -  0.d+00 
deta  -  delu 
dchi  »  0.d+00 
dpa  -  0 . d+00 
elseif(n.eq.3)  then 
dthe  «  0.d+00 
deta  -  0.d+00 
dch i  -delu 
dpa  -  0 . d+00 
elsei f(n.eq.4)  then 
dthe  -  0.d+00 
deta  -  0.d+00 
dch i  -  0 . d+00 
dpa  -delu 
endi  f 
C 

COHMEHT:  The  folloeing  updates  the  angles  in  DELTR  for  PU  eatrix 
c  determination.  The  suffix  "p*  (plus)  means  that  the  del 

c  I s  added  ehen  app I i cab  I e . 

c 


thetap  -  theta  +  dthe 


comment:  This  following  section  sets  up  the  emus  delu  s. 
c  It  updates  the  angles  in  DELTR  for  PU  determination, 

c  The  suffix  ■  means  that  the  delu  is  subtracted  mhen 

c  applicable.  HOTE:  The  dthe,deta, . . .  are  negative  nom 


comment:  This  2nd  call  to  DELTA  calculates  the  F(minus)  values. 


cal  I  DELTR(a, i jmu, thetam^etam^chim^pam, tpi ,sac,tp,dam) 


comment:  The  fol losing  sets  up  the  approximations  of  each  partial 
c  Print  out  these  (dap  -  dam)  values: 

c  The  fol  losing  sets  up  flatrix  PU  [1x1].  Note  the 

c  denominator. This  is  2  x  delu  because  of  the  method  of 

c  approx i mat i on . 

c 

dif  -  dap  -  dam 

if(part.eq.l)  C(n)  -  di f/(2.d+00*delu) 
i f(part .eq.2)  pu(s;n)  -  di f/(2.d+00*delu) 
c 

comment:  start  forming  the  elements  of  the  PU  matrices, 
c  if  k  -  1,  then  PUB  is  being  formed, 

c  if  k  ■  2,  then  PUB  is  being  formed. 


100  cont i nue 
return 


c _ 

SUBROUT  I  HE  DELTR(a,  i ,  mu,  theta,  et  a,  ch  i ,  pa,  tp  i ,  sac, tp, da) 
c 

comment:  This  subroutine  calculates  the  orbital  parameter 
c  perturbation  due  to  changes  in  orientation  angles 

c  (theta, eta, ch i , pa)  and  changes  in  orbital  parameters, 

c 

DOUBLE  PRECISION  01 ,02,03, Bt ,B2,B3,B-t,B5,B6,B?,B8 
DOUBLE  PRECISION  theta, eta, chi ,pa,tp,tpi ,a, i , mu, sac 
DOUBLE  PRECISION  da, dal ,da2,da3,da1,da5,da6 


C. 


B1 

■  dsin(theta)*dcos(chi ) 

B2 

■  dsin(theta)*dsin(chi )*dcos(eta) 

B3 

■  dcos(theta)*dsin(chi )*dsin(eta) 

B1 

*  dsin(theta)*dsin(chi ) 

B5 

■  dsin(theta)*dcos(eta) 

B6 

■  dcos(theta)*dcos(chi )*dsin(eta) 

B7 

-  dsin(theta)*dsin(eta) 

B8 

n 

■  dcos(theta)*dcos(eta) 

01 

-  dsin(theta)*dcos(chi )*dcos( i ) 

D2 

-  dsin(theta)*dsin(chi )*dcos(eta)*dcos( i ) 

$ 

♦  dsin(theta)*dsin(eta)*dsin( i ) 

03 

■  dcos(theta)*dsin(chi )*dsin(eta)*dcos( i ) 

$  -  dcos(theta)*dcos(eta)*dsin( i ) 


c. _ _ _ 

comment:  Calculate  the  fo Homing  factors  that  enter  into  the  0a 
c  equat i on . 


dal  -  01 **2 . d+00/1 . d+00*( (B2  ♦  3.d+00*B1)*dcos(pa) 

$  -  (3.0+00*B1  ♦  B5)*dsin(pa)) 

da2  -  -01*02/2. D+00*((B2-B1)*dsin(pa)  ♦  (B5  -  B1 )*dcos(pa)) 
da3  -  01 *03*2 . d+00*(B6*dcos(pa)  -  B3*dsin(pa)) 
da4  -  02*D3*2.d+00*(B3*dcos(pa)  +  B6*dsin(pa)) 
da5  -  D2**2 . d+00/1 . d+00*( ( 3 . d+00*B2  ♦  B4)*dcos(pa) 

$  -  (B1+3.d+00*85)*dsin(pa)) 

da6  -  03**2 . d+00*( (B2+B4)*dcos(pa)  -  (B1  ♦  B5)*dsin(pa)) 
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couent:  Calculte  the  change  in  Seai-aajor  axis, 
c 

da  *  a**1 .5d+00*sac*tp/iu*(da1+da2+da3+da4+da5+da6) 
c 

■rite(18,581 )  da 

581  fonat(/5x,‘Da  -  ”,  f20. 10) 


return 


SUBROUTINE  ECHO(ax,  ix,tp,etax,chix,thetax,pax) 

THIS  SUBROUTINE  ECHOES  IHPUT  TO  SCREEN  FOR  UERI FI  CATION. 

DOUBLE  PRECISION  ax, ix,tp,etax,chix,thetax,pax 
•rite(*,500) 

■ri te(1 9, 500) 

•rite(*,501)ax, ix,tp 
•rite(19,501)ax, ix,tp 
erite(*,502)etax,chix,thetax,pax 
•r i te( 1 9, 502 )et ax, ch i x, thetax,  pox 

_  FORNRTS  _ 


500  FORMAT (1  OX, 'YOUR  SYSTEM  IS  RS  FOLLOWS :*) 

501  F0RMRT(20X, ’ORBITAL  PARAMETERS:  a  (DU)  -’,017.10, 

S/20X, '  i  (DEG)  -’,017.10, 

$/20X, '  TP  (TU)  -’,017.10) 

502  F0RMRT(20X,'H  ORIENTATION:  eta  (DEG)  -',017.10, 

$/20X, '  chi  (DEG)  — ’ ,017. 10, 

$/20X, '  CONING  ANGLE:  theta  (DEG)  -',017.10, 

$/20X , '  PHASE  ANGLE  aa  (DEG)  -’ .017.10) 


SUBROUTINE  INU4X4(PC,PCI  ^et ) 

C 

COnUENT:  This  subroutine  calculates  the  inuerse  of  the  [4x1]  PC 
C  eatrix  using  CORNER'S  RULE.  Primitive,  but  it  eorks! 

C _ 

double  precision  pc(4,4),pci(4,4),det 
double  precision  dpi  1  Jdp12,dpl3JdpHJdp21  ,dp22,dp23,dp21 
double  precision  dp31 ,dp32,dp33,dp34,dp41 ,dp12,dp43,dp41 
double  precision  d1122id1123ld1124ld1221id1223ld1221 
double  precision  dl 321  ,d1322,d1321,d1421  ,d1422,d1423 
double  precision  d2112^d2113id21MJd2211>d2213^d22M 
double  precision  d231 1  ,d2312,d2314,d241 1  ,d2412,d2113 
double  precision  d31 12,d3113jd31 14,d3211,d3213,d3214 
double  precision  d331 1  ,d3312,d3314,d341 1  ,d3412,d3413 
double  precision  dll  12,d41 13, dll 11,d4211  ,d121 3,d4211 
double  precision  d431 1  ,d431 2 ,d4314,d441  1  fd1412,d4413 


COIfflENT: 

c 

f! 

Deteraine  the  eleaents  of  eoch  roa  of  the  deterainant 
and  cofactor  Matrix.  Data  entered  by  colunn. 

COUtlEHT : 

r 

1st  roa  eleaents: 

n 

dl 122  -  pc(3,3)*pc(1,4)  -  pc(4,3)*pc(3,4) 
dl 1 23  *  pc(3,2)*pc(4,1)  -  pc(1,2)*pc(3,4) 
dl 1 24  ■  pc(3,2)*pc(4,3)  -  pc(4#2)*pc(3,3) 
dpi  1  »  pc(2,2)*d1 122  -  pc(2,3)*d1123  ♦  pc(2,1)*d1124 

c _ 

dl 221  -  dl 122 

dl 223  -  pc(3, 1 )*pc(4,1)  -  pc(1#1)*pc(3,1) 
d1224  -  pc(3,1)*pc(4,3)  -  pc(4, 1)*pc(3,3) 
dpi 2  -  pc(2,1)*d1221  -  pc(2,3)*d1223  ♦  pc(2,1)*d1221 

dl 321  -  dl 123 
dl 322  -  d 1 223 

dl 324  -  pc(3,1)*pc(1,2)  -  pc(4, 1 )*pc(3,2) 

dpi 3  -  pc(2,1)*d1321  -  pc(2,2)*d1322  ♦  pc(2,1)*d1324 
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dH21  -  dll 24 
dH22  -  d1224 
dl 423  -  dl 324 

dpi 4  -  pc(2,1)*d1421  -  pc(2,2)*d1422  ♦  pc(2#3)*d1423 

*  C _ 

COHHENT:  2nd  ro«  elements; 

C _ 

d21 1 2  ■  pc(3,3)*pc(4,4)  -  pc(4,3)*pc(3,4) 

*  d2113  -  pc(3,2)*pc(4,4)  -  pc(4,2)*pc(3,4) 

d2114  *  pc(3,2)*pc(4,3)  -  pc(4,2)*pc(3,3) 

dp21  -  pc(1 ,2)*d21 12  -  pc(1 ,3)*d21 13  ♦  pc(1 ,4)*d21 14 

C _ 

d2211  -  d21 12 

d2213  -  pc(3j 1 )*pc(4,4)  -  pc(4, 1 )*pc(3,4) 
d2214  -  pc(3J)*pc<4,3)  -  pc(4,1)*pc(3,3) 
dp22  -  pc(1,1)*d2211  -  pc(1 ,3)*d2213  ♦  pc( 1 , 4)*d221 4 
C _ _ _ 

*  d23J 1  -  d2J13 

d2312  -  d2213 

d2314  -  pc(3, 1 )*pc(4,2)  -  pc(4 ,1 )*pc(3,2) 

dp23-  pc(1,1)*d2311  -  pc(1,2)*d2312  ♦  pc(1 ,4)*d2314 

*  C - 

d2411  -  d2114 
d2412  -  d2214 
d2413  -  d231 4 

dp24  -  pc(1J)*d2411  -  pc(1 ,2)*d2412  ♦  pc(1,3)*d2413 

*  C _ 

COHUEHT :  3rd  rov  elements; 

C _ 

d31 12  -  pc(2,3)*pc<4,4)  -  pc(4,3)*pc(2,4) 

%  d3113  -  pc(2,2)*pc(4,4)  -  pc(4,2)*pc(2,4) 

d31 14  -  pc(2,2)*pc(4,3)  -  pc(4#2)*pc(2,3) 

dp31  -  pc(1,2)*d3112  -  pc(1 ,3)*d3113  ♦  pc(1 ,4)*d31 14 


d321 1  -  d31 1 2 

d3213  -  pc(2,1)*pc(1,4)  -  pc(4l1)*pc(2#4) 
d3214  •  pc(2,1)*pc(4,3)  -  pc(4,1)*pc(2,3) 
dp32  -  pc(1J)*d3211  -  pc(1,3)*d3213  ♦  pc(M)*d32H 

d331 1  -  d31 1 3 
d3312  -  d321 3 

d3314  -  pc(2,1)*pc(4,2)  -  pc(4, l)*pc(2,2) 

dp33  -  pc(1 , 1  )*d331 1  -  pc(1 ,2)*d3312  ♦  pc( 1,1)*d3311 


d341 1  -  d3114 
d3412  -  d32M 
d3413  -  d33 14 
dp34  -  pc(1,1)*d3411 
r 


COntlEHTS:  1th  ro»  elements: 


-  pc(1#2)*d3412+  pc(1,3)*d3413 


dll  12  -  pc(2,3)*pc(3,4)  -  pc(3,3)*pc(2,1) 
dll  13  -  pc(2,2)*pc<3,4)  -  pc(3,2)*pc(2#1) 
dllll  -  pc(2j2)*pc(3,3)  -  pc(3,2)*pc(2,3) 
dpll  -  pc(1,2)*d1112  -  pc(l,3)*d4H3  ♦  pc(M)*d41M 


d1211  -  dll  12 

d4213  -  pc(2,1)*pc(3,4)  -  pc(3# 1 )*pc(2,1) 
d1214  -  pc(2,1)*pc(3#3)  -  pc(3, 1 )*pc(2,3) 
dp42  -  pc(1,1)*d1211  -  pc(l,3)*d4213  ♦  pc(1 ,4)*d4214 


d431 1  -  dll  13 
d4312  -  d4213 

d4314  »  pc(2,1)*pc(3,2)  -  pc(3, 1 )*pc(2,2) 

dp43  -  pc(1  J)*d1311  -  pc(1,2)*d4312  ♦  pc(1,1)*d4314 


d441 1  -  dllll 
dll  12  -  d42H 
d4413  -  d1314 

dpll  -  pc(1l1)*d4411  -  pc(1,2)*d4412  ♦  pc(1 ,3)*d4413 


tS' 


■*.  f ,  '  *  **. 
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C0W1ENT:  Calc  the  determinant  of  [PC]. 


det  ■  pc(1,1)*dp11  -  pc(1,2)*dp12  ♦  pc(1,3)*dp13  - 
$  PC(1,4)*dp14 


C0W1ENT:  Flag  the  case  ehen  a  singular  eatrix  exist., i.e.,  det*  0.0 
C 

i f (det .eq.0.d+00)  then 
print*, 'DANGER!  DANGER !  Det  [PC]  -  0.0' 
print*, 'Prograe  is  stopped  in  Subroutine  INU4X4:no  output' 
go  to  208 
end 


CONNENT:  Calculate  the  Inverse  [Pci]  of  [Pc]. 


1,2 

1,3 

M 

2,1 

pci (2,2) 
pci (2,3) 
pci (2,4) 
pci (3,1) 
pci (3,2) 
pci(3,3) 
pci (3, 4) 
pci (4,1) 
pci(4,2) 
pci (4,3) 
pci(4,4) 

cont i nue 
return 


dpi  1 /det 

21, 

31, 

41 

12, 

dp22/det 
-dp32/det 
dp42/det 
dpi  3/de t 
-dp23/det 
dp33/det 
-dp43/det 
-dp14/det 
dp24/det 
-dp34/det 
dp44/det 


i 


Attttfinda  I 

nptimgf1if>n  Hal  &ge  Output 

A  test  case  vu  run  to  determine  the  required  changes  in  an  initial  set  of  control 
parameter  necessary  to  reach  an  stationary  value  for  the  Performance  Index  J.  The 
specific  input  data  with  the  resulting  output  are  given  below.  For  purpose  of 
demonstration,  two  iterations  are  shown  here.  An  option  does  exist  for  n  iterations  or 
until  a  certain  predetermined  convergence  value  is  reached.  Convergence  is  set  at  .001 
radians  or  .037  degrees  and  is  reached  in  two  iterations. 

Case  Tested.  The  case  (depicted  by  Figure  2 AD)  is  chosen  as  the  test  for  the 
search  scheme  performance  demostration.  The  following  are  the  input  data: 

a  -  2.0  du 

i  -  *5  0  * 

8  -  30* 

T)  -  0.0* 

<  -  0.0- 

PA  -  -90* 

From  the  Figure  2.4C.  the  local  maximum  appears  to  be  about  33*  This  is  what  is 
expected  from  the  search  scheme.  As  is  shown  later  in  the  next  few  pages,  the 
perspective  provided  a  very  good  estimate  of  where  the  actual  maximum  is  located  with 
the  initial  conditions  given  to  within  .3  degrees.  This  is  not  always  the  ca at  since 
pinpointing  the  aiaximum  from  the  figures  is  limited  to  the  amount  of  divisions  in  the 

two  variables  8  and  PA.  However,  increments  of  3  degrees  is  sufficient  to  make  a  good 
initial  guess  with.  This  is  strongly  exemplified  by  the  test  case  which  follows: 
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SYSTEM  IS  AS  FOLLOUS: 


ORBITAL  PARAMETERS:  a  (DU) 

i  (DEG) 

tp  (TU) 

ORIENTATION:  eta  (DEG) 

chi  (DEG) 

CONING  ANGLE:  theta  (DEG) 

PHASE  ANGLE:  pa  (DEG) 


.2000000Q00e+01 
.4500000000e+02 
. 17771 531 75e+02 
. OOOOOOOOOOe+OO 
. OOOOOOOOOOe+OO 
. 5000000000e+02 
-.9000000000e+02 


Uector  C  [1x4] 

.0000180033  .0000000000  .0000000000 


.0000000000 


For  k  ■  1,  Matrix  PUA  [4x4] 

. 65971 25567e-05  .OOOOOOOOOOe+OO  .OOOOOOOOOOe+OO  .OOOOOOOOOOe+OO 
. 1 836358803e-04  . 271 9543322e-05  - . 21 56448749e-05  - . 3691 1 1 4090e-05 
.l795098150e-04  -. 21 40999434e-05  -.381 1 10577 9e- 05  ,6321384922e-05 
. 1 798084043e-04  - . 36971 25027e-05  . 6322957922e-05  - . 6323023559e-05 

For  k  ■  2,  Matrix  PUB  [4x4] 

.3087931771e-04  .OOOOOOOOOOe+OO  .OOOOOOOOOOe+OO  .OOOOOOOOOOe+OO 
. 1 836358803e-04  -.2719543322e-05  .2156448749e-05  .3691 1 14090e-05 
.1795098150e-04  .2140999434e-05  .381 1 1 05779e-05  -.6321384922e-05 
. 1 798084043e-04  . 36971 25027e-05  - . 6322957922e-05  . 6323023559e-05 


Matrix  DIFF  [4x4] 

-.242821 921 4e-04  .OOOOOOOOOOe+OO  .OOOOOOOOOOe+OO  .OOOOOOOOOOe+OO 
.OOOOOOOOOOe+OO  . 5439086644e-05  -,4312897499e-05  -.73822281 81 e-05 
. OOOOOOOOOOe+OO  - . 4281 998867e-05  - . 762221 1 558e-05  . 1 264276984e-04 
. OOOOOOOOOOe+OO  - . 7394250055e-05  . 1 264591 584e-04  - . 1 26460471 2e-04 

Matrix  PC  [4x4] 

.2428219214e-03  .OOOOOOOOOOe+OO  .OOOOOOOOOOe+OO  .OOOOOOOOOOe+OO 
. OOOOOOOOOOe+OO  - . 5439086644e-04  . 431 289?499e-04  . 73822281 81 e-04 
.OOOOOOOOOOe+OO  .4281998867e-04  .762221 1558e-04  - . 1 264276984e-03 
.OOOOOOOOOOe+OO  . 7394250055e-04  -.1 26459 1584e-03  .1 2646047 12e-03 


» 


Deters  inant  -  -.000000000000000268858 1996521702033557475 


Hatrix  PCI  [4x4] 

.41 18244325e+04  . 0000000000e+00  . 0000000000e+00  . 0000000000e+00 

>  . 0000000000e+00  . 5734031 325e+04  . 133573851 4e+05  . 1000663821 e+05 

:  . 0000000000e+00  . 13333724 13e+05  . 1 1 142191 14e+05  . 335563861 6e+04 

. 0000000000e+00  . 9980849403e+04  .3331 900 100e+04  .541 22476 I6e*04 

Uector  C  [1x4] 

>  . 1800333992e-04  . 0000000000e+00  . 0000000000e+00  . OOOOOOOOOOe+Ol 

Uector  DEL  (rads)  [1x4] 

.7414215246e-01  . 0000000000e+00  . 0000000000e+00  . 0000000000e+00 

Uector  DEL  (degs)  [1x4] 

.4248032420e+01  . 0000000000e+00  . 0000000000e+00  .0000000000e+00 


This  is  the  end  of  the  INITIAL  LOOP. 
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ITERATION  -  1 


SYSTEfl  IS  AS  FOLLOUS: 


ORBITAL  PflRflflETERS : 

a  (DU) 

s 

.2000000000e+01 

i  (DEG) 

m 

. 1500000000e+02 

TP  (TU) 

m 

. 17771 531 75c+02 

ORIENTATION: 

eta 

(OEG) 

m 

. 0000000000e+00 

chi 

(OEG) 

m 

. 0000000000e+00 

CONING  ANGLE: 

theta  (DEG) 

m 

. 5121803212e+02 

PHASE  ANGLE 

PR 

(DEG) 

m 

- . 9000000000e+02 

Oector  C  [1x1] 

.0000016595  .0000000000  .0000000000  .0000000000 


For  k  *  1,  flatrix  PUR  [1x1] 
-.7317116303e-05  . 0000000000e+00  . 0000000000e+00 

. 20375 13771e-05  .3816550277e-05  - . 2113709575e-06 
. 161 655571 6c-05  -.227391 1 101e-06  -.395 15532 IOc-05 
. 1657156877e-05  -.3236116172e-05  . 6358261017e-05 


. 0000000000e+00 
- . 3230390882e-05 
. 63563937 11e-05 
-.6358312092c-05 


For  k  -  2.  flatrix  PUB  [1x1] 


.1230272911e-01 
.  203751 3771e-05 
.1616555716e-05 
.  1 657156877e-05 


. 0000000000e+00 
-.3816550277e-05 
.2273911101e-06 
.3236116172e-05 


.0000000000e+00 
. 2113709575e-06 
.3951 55321 8e-05 
- . 6358261017e-05 


. 0000000000e+00 
.32303908826-05 
-. 635639371 1e-05 
. 6358312092e-05 


flatrix  DIFF  [1x1] 


- . 1 96201 7511 e-01  . 0000000000e+00 
. 0000000000e+00  . 7633 1 00551e-05 
. 0000000000e+00  -.1517822207e-06 
. 0000000000e+00  - . 6172292913e-05 


.  0000000000e+00  . 0000000000e+00 
.18871191506-06  -. 616078 1761e-05 
.  79031 06136e-05  . 1271278713e-01 
.1271 652809e-01  -.1271 6681 1 8e-01 


.  196201 7511 e-03 
. 0000000000e+00 
. 0000000000e+00 
. 0000000000e+00 


flatrix  PC 
. 0000000000e+00 
-.7633100551e-01 
.15178222076-05 
. 6172292913c-01 


[1x1] 

. 0000000000e+00 
.18871191506-05 
.7903106136e-01 
.12716528096-03 


.00000000006+00 
. 6160781 761e-01 
-. 1271278713e-03 
. 1271 66811 8e-03 


s 


•  Determinant  -  .0000000000000000109809814800276459131665 

Matrix  PCI  [4x4] 

. 5096794392e+04  . 0000000000e+00  . 0000000000e+00  . 0000000000e+00 

.QOOOOOOOOOe+OO  - . 1 092792724e+06  -.1579012964e+06  - . 1 023329720e+06 

•  . QOOOOOOOOOe+OO  - . 1 573479727e+06  -.2481494526e+06  -. 1681319044e+06 
. 0000000000e+00  -.1017271 824e+06  - . 1 677808467e+06  -.1081 826908e+06 

Uector  C  [1x4] 

•  . 1 659530858e-05  . 0000000000e+00  . 0000000000e+00  . 0000000000e+00 

Uector  DEL  (rads)  [1x4] 

. 8458287572e-02  .0000000000e+00  . 0000000000e+00  . 0000000000e+00 

**  Uector  DEL  (degs)  [1x1] 

. 484621 1798e+00  . OOOOOOOOOOe+OO  . 0000000000e+00  .OOOOOOOOOOe+OO 


ITERATION  -  2 


SVSTEB  IS  AS  FOLLOUS: 


ORBITAL  PRRRBETERS: 

a  (DU) 

s 

, 2000000000e+01 

i  (DEG) 

m 

. 1500000000e+02 

TP  (TU) 

- 

. 17771 531 75e+02 

ORIENTATION: 

eta 

(DEG) 

m 

. 0000000000e+00 

chi 

(DEG) 

m 

. 0000000000e+00 

CONING  ANGLE: 

theta  (DEG) 

m 

. 5173265660e+02 

PHASE  ANGLE 

PA 

(DEG) 

3 

- . 9000000000e+02 

Uector  C  [1x1] 

.0000000211  .0000000000  .0000000000  .0000000000 


Fop  k  -  1,  Batrix  PUA  [1x1] 

- . 8661821950e-05  . 0000000000e+00  .0000000000e+00  .0000000000e+00 
,1001718609e-06  . 39111971 62e-05  1718988786e-07  -.31 7373288 1e-05 

- . 2087678209e-07  - . 3857831 379c-1 0  - . 3965810001 e-05  . 635670691 7e-05 

.2105726398e-07  -.31 7915167 I c-05  .635861 1121c-05  - . 6358690891e-05 


For  k  -  2.  tlatrix  PUB  [1x1] 


. 1 03899351 8e-01  . 0000000000e+00 
.1001718609e-06  -. 3911197 162e-05 
- . 2087678209c-07  . 3857831 379e- 1 0 
.2105726398e-07  .31 79151671 e-05 


. 0000000000e+00  . 0000000000e+00 

. 1718988786e-07  .31 73732881 e-05 

.3965810001 e-05  -. 63567069 17e-05 
- . 635861 1 121e-05  . 6358690891e-05 


Batrix  DIFF  [1x1] 

- . 1 90517601 3e-01  . 0000000000e+00  . 0000000000e+00 

. 0000000000e+00  . 7888991325e-05  -. 3137977571 e-07 
. 0000000000e+00  -. 771 5662756c- 1 0  -.7931680002e-05 
. 0000000000e+00  - . 6358909313e-05  .1271 722285e-01 


. 0000000000e+00 
- . 6317165762e-05 
. 127 1311 383e-01 
-. 1271 7381 79e-01 


Batrix  PC 

. 1 90517601 3e-03  . 0000000000e+00 
. 0000000000e+00  - . 7888991325e-01 
. 0000000000e+00  ,7715662758e-09 


[1x1] 

. 0000000000e+00 
. 3137977571 e-06 
,7931680002e-01 


. 0000000000e+00  . 6358909313e-01  -.1271 722285e-03 


.0000000000e+00 
. 6317165762e-01 
127 1311 383e-03 
. 1271 7381 79e-03 


Oetera i nant  -  . 000000000000000029876 ! 55962369095986375 1 


flat r i x  PCI  [4x1] 

. 5248032476e+04  . OOOOOOOOOOe+OO  . OOOOOOOOOOe+OO  .OOOOOOOOOOe+OO 
. OOOOOOOOOOe+OO  - . 3878353975e+05  -.51 7627531 7e+05  -.3238906519e+05 
. OOOOOOOOOOe+OO  -.5156189311e+05  -.89?3109370e+05  -. 6396766 162e+05 
. OOOOOOOOOOe+OO  -.3216881292e+05  -.6384770652e+05  -.3990851 93 1e+05 


Uector  C  [1x4] 

.21 08361 300e-07  .OOOOOOOOOOe+OO  .OOOOOOOOOOe+OO  .OOOOOOOOOOe+OO 


Uector  DEL  (rads)  [1x1] 

. 1 1 06474857e-03  .OOOOOOOOOOe+OO  .OOOOOOOOOOe+OO  .OOOOOOOOOOe+OO 


Uector  DEL  (degs)  [1x4] 

. 6339633946e-02  .OOOOOOOOOOe+OO  .OOOOOOOOOOe+OO  .OOOOOOOOOOe+OO 


***********************************************s*»**«******«******** 


By  jove,  I  think  you  got  it.  Solar  Sailor!!!!! 

******************************************************************** 


FIHBL  SYSTEH  IS  RS  FOLLOUS: 

ORB  URL  PRRflflETERS :  a  (DU) 

i  (DEG) 


OR  I EHTRT I  OH : 


CORING  RHGLE : 
PHRSE  RNGLE 


tp  (TU)  - 
eta  (DEG)  - 
chi  (DEG)  - 
theta  (DEG)  * 
pa  (DEG)  - 


.2000000000e+01 
. 1500000000e+02 
. 17771 531 75e+02 
.OOOOOOOOOOe+OO 
.OOOOOOOOOOe+OO 
. 5473899623e+02 
- . 9000000000e+02 


This  is  the  end  of  the  SECOND  ITERATION. 
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REtlRRKS:  This  test  case  shoaed  that  ahen  the  search  point 
•as  initiated  at  theta  -  50° ,  the  search  algorithm  converged  at  a 
local  sax i bub  at  theta  *  54.73°.  Froa  Figure  4.2C,  this  can  be 
visual Ig  verified.  The  close  agreement  between  the  perspectives  and 
the  computed  box i bub  gives  strong  credence  in  the  aethod  of  search 
and  the  utility  of  the  perspectives  in  providing  a  good  initial 
guess.  This  last  control  vector  set  aas  inputted  into  the  DELTA  A 
function  to  verify  the  finding  aithin  the  theta  equal  54.7®.  The 
results  are  tabulated  beloa  for  a  fix  inclination  at  45°: 

THETA  DELTA  A  (x  10"4  DU) 


54.0 

1 .272  129 

54.1 

1.272  170 

54.2 

1.272  205 

54.3 

1.272  233 

54.4 

1.272  256 

54.5 

1.272  273 

54.6 

1.272  284 

— ►  54.7 

1.272  289 

54.8 

1.272  288 

54.9 

1.272  281 

55.0 

1.272  269 

The  search  scheae  is  aore  accurate  in  pin-pointing  the  exact  aaxieua 
coordinate.  The  coordinate  is  essentially  the  control  vector  set 
that  aould  optiaize  the  change  in  seai-aajor  axis.  Note,  hoaever, 
that  other  aaxiaa  do  exist  and  can  easily  be  located  aith  the  aid  of 
the  correspond  I ng  perspective  and  /or  aith  the  search  a  Igor  it ha. 
C are  aust  be  excercised  ahen  doing  so;  the  perspectives  give  the 
angular  aoaentua  orientation  angles  as  zero  degrees.  The  search 
algoritha  ail  I  find  aaxiaa  aithout  holding  these  angles  fixed.  All 
variables  are  Increaented.  Plotting  the  perspective  for  a 
particular  output  ail  I  provide  an  extra  diaension  in  locateing  and 
verifying  the  coaputed  aaxiaa. 
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